newmatrm.cpp
Go to the documentation of this file.00001
00002
00003
00004
00005 #ifndef WANT_MATH
00006 #define WANT_MATH
00007 #endif
00008
00009 #include "newmat.h"
00010 #include "newmatrm.h"
00011
00012 #ifdef use_namespace
00013 namespace NEWMAT {
00014 #endif
00015
00016 #ifdef DO_REPORT
00017 #define REPORT { static ExeCounter ExeCount(__LINE__,12); ++ExeCount; }
00018 #else
00019 #define REPORT {}
00020 #endif
00021
00022
00023
00024
00025
00026 void RectMatrixRow::Reset (const Matrix& M, int row, int skip, int length)
00027 {
00028 REPORT
00029 RectMatrixRowCol::Reset
00030 ( M.Store()+row*M.Ncols()+skip, length, 1, M.Ncols() );
00031 }
00032
00033 void RectMatrixRow::Reset (const Matrix& M, int row)
00034 {
00035 REPORT
00036 RectMatrixRowCol::Reset( M.Store()+row*M.Ncols(), M.Ncols(), 1, M.Ncols() );
00037 }
00038
00039 void RectMatrixCol::Reset (const Matrix& M, int skip, int col, int length)
00040 {
00041 REPORT
00042 RectMatrixRowCol::Reset
00043 ( M.Store()+col+skip*M.Ncols(), length, M.Ncols(), 1 );
00044 }
00045
00046 void RectMatrixCol::Reset (const Matrix& M, int col)
00047 {
00048 REPORT
00049 RectMatrixRowCol::Reset( M.Store()+col, M.Nrows(), M.Ncols(), 1 );
00050 }
00051
00052
00053 Real RectMatrixRowCol::SumSquare() const
00054 {
00055 REPORT
00056 long_Real sum = 0.0; int i = n; Real* s = store; int d = spacing;
00057
00058 if (i) for(;;)
00059 { sum += (long_Real)*s * *s; if (!(--i)) break; s += d; }
00060 return (Real)sum;
00061 }
00062
00063 Real RectMatrixRowCol::operator*(const RectMatrixRowCol& rmrc) const
00064 {
00065 REPORT
00066 long_Real sum = 0.0; int i = n;
00067 Real* s = store; int d = spacing;
00068 Real* s1 = rmrc.store; int d1 = rmrc.spacing;
00069 if (i!=rmrc.n)
00070 {
00071 Tracer tr("newmatrm");
00072 Throw(InternalException("Dimensions differ in *"));
00073 }
00074
00075 if (i) for(;;)
00076 { sum += (long_Real)*s * *s1; if (!(--i)) break; s += d; s1 += d1; }
00077 return (Real)sum;
00078 }
00079
00080 void RectMatrixRowCol::AddScaled(const RectMatrixRowCol& rmrc, Real r)
00081 {
00082 REPORT
00083 int i = n; Real* s = store; int d = spacing;
00084 Real* s1 = rmrc.store; int d1 = rmrc.spacing;
00085 if (i!=rmrc.n)
00086 {
00087 Tracer tr("newmatrm");
00088 Throw(InternalException("Dimensions differ in AddScaled"));
00089 }
00090
00091 if (i) for (;;)
00092 { *s += *s1 * r; if (!(--i)) break; s += d; s1 += d1; }
00093 }
00094
00095 void RectMatrixRowCol::Divide(const RectMatrixRowCol& rmrc, Real r)
00096 {
00097 REPORT
00098 int i = n; Real* s = store; int d = spacing;
00099 Real* s1 = rmrc.store; int d1 = rmrc.spacing;
00100 if (i!=rmrc.n)
00101 {
00102 Tracer tr("newmatrm");
00103 Throw(InternalException("Dimensions differ in Divide"));
00104 }
00105
00106 if (i) for (;;) { *s = *s1 / r; if (!(--i)) break; s += d; s1 += d1; }
00107 }
00108
00109 void RectMatrixRowCol::Divide(Real r)
00110 {
00111 REPORT
00112 int i = n; Real* s = store; int d = spacing;
00113
00114 if (i) for (;;) { *s /= r; if (!(--i)) break; s += d; }
00115 }
00116
00117 void RectMatrixRowCol::Negate()
00118 {
00119 REPORT
00120 int i = n; Real* s = store; int d = spacing;
00121
00122 if (i) for (;;) { *s = - *s; if (!(--i)) break; s += d; }
00123 }
00124
00125 void RectMatrixRowCol::Zero()
00126 {
00127 REPORT
00128 int i = n; Real* s = store; int d = spacing;
00129
00130 if (i) for (;;) { *s = 0.0; if (!(--i)) break; s += d; }
00131 }
00132
00133 void ComplexScale(RectMatrixCol& U, RectMatrixCol& V, Real x, Real y)
00134 {
00135 REPORT
00136 int n = U.n;
00137 if (n != V.n)
00138 {
00139 Tracer tr("newmatrm");
00140 Throw(InternalException("Dimensions differ in ComplexScale"));
00141 }
00142 Real* u = U.store; Real* v = V.store;
00143 int su = U.spacing; int sv = V.spacing;
00144
00145
00146
00147
00148
00149 if (n) for (;;)
00150 {
00151 Real z = *u * x - *v * y; *v = *u * y + *v * x; *u = z;
00152 if (!(--n)) break;
00153 u += su; v += sv;
00154 }
00155 }
00156
00157 void Rotate(RectMatrixCol& U, RectMatrixCol& V, Real tau, Real s)
00158 {
00159 REPORT
00160
00161 int n = U.n;
00162 if (n != V.n)
00163 {
00164 Tracer tr("newmatrm");
00165 Throw(InternalException("Dimensions differ in Rotate"));
00166 }
00167 Real* u = U.store; Real* v = V.store;
00168 int su = U.spacing; int sv = V.spacing;
00169
00170
00171
00172
00173
00174
00175 if (n) for(;;)
00176 {
00177 Real zu = *u; Real zv = *v;
00178 *u -= s * (zv + zu * tau); *v += s * (zu - zv * tau);
00179 if (!(--n)) break;
00180 u += su; v += sv;
00181 }
00182 }
00183
00184
00185
00186
00187 Real pythag(Real f, Real g, Real& c, Real& s)
00188
00189
00190
00191 {
00192 if (f==0 && g==0) { c=1.0; s=0.0; return 0.0; }
00193 Real af = f>=0 ? f : -f;
00194 Real ag = g>=0 ? g : -g;
00195 if (ag<af)
00196 {
00197 REPORT
00198 Real h = g/f; Real sq = sqrt(1.0+h*h);
00199 if (f<0) sq = -sq;
00200 c = 1.0/sq; s = h/sq; return sq*f;
00201 }
00202 else
00203 {
00204 REPORT
00205 Real h = f/g; Real sq = sqrt(1.0+h*h);
00206 if (g<0) sq = -sq;
00207 s = 1.0/sq; c = h/sq; return sq*g;
00208 }
00209 }
00210
00211
00212
00213
00214 #ifdef use_namespace
00215 }
00216 #endif
00217
00218