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