Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

newmatrm.cpp

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

newmat11b
Generated Mon May 9 04:54:18 2016 by Doxygen 1.6.3