newmatrm.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 // Copyright (C) 1991,2,3,4: R B Davies
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 // operations on rectangular matrices
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    // while (i--) { sum += (long_Real)*s * *s; s += d; }
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    // while (i--) { sum += (long_Real)*s * *s1; s += d; s1 += d1; }
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    // while (i--) { *s += *s1 * r; s += d; s1 += d1; }
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    // while (i--) { *s = *s1 / r; s += d; s1 += d1; }
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    // while (i--) { *s /= r; s += d; }
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    // while (i--) { *s = - *s; s += d; }
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    // while (i--) { *s = 0.0; s += d; }
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    //while (n--)
00147    //{
00148    //   Real z = *u * x - *v * y;  *v =  *u * y + *v * x;  *u = z;
00149    //   u += su;  v += sv;
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    //  (U, V) = (U, V) * (c, s)  where  tau = s/(1+c), c^2 + s^2 = 1
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    //while (n--)
00172    //{
00173    //   Real zu = *u; Real zv = *v;
00174    //   *u -= s * (zv + zu * tau); *v += s * (zu - zv * tau);
00175    //   u += su;  v += sv;
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 // misc procedures for numerical things
00188 
00189 Real pythag(Real f, Real g, Real& c, Real& s)
00190 // return z=sqrt(f*f+g*g), c=f/z, s=g/z
00191 // set c=1,s=0 if z==0
00192 // avoid floating point overflow or divide by zero
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;           // make return value non-negative
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 


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:13