newmatrm.cpp
Go to the documentation of this file.
1 
6 
7 // Copyright (C) 1991,2,3,4: R B Davies
8 
9 #define WANT_MATH
10 
11 #include "newmat.h"
12 #include "newmatrm.h"
13 
14 #ifdef use_namespace
15 namespace NEWMAT {
16 #endif
17 
18 #ifdef DO_REPORT
19 #define REPORT { static ExeCounter ExeCount(__LINE__,12); ++ExeCount; }
20 #else
21 #define REPORT {}
22 #endif
23 
24 
25 // operations on rectangular matrices
26 
27 
28 void RectMatrixRow::Reset (const Matrix& M, int row, int skip, int length)
29 {
30  REPORT
32  ( M.Store()+row*M.Ncols()+skip, length, 1, M.Ncols() );
33 }
34 
35 void RectMatrixRow::Reset (const Matrix& M, int row)
36 {
37  REPORT
38  RectMatrixRowCol::Reset( M.Store()+row*M.Ncols(), M.Ncols(), 1, M.Ncols() );
39 }
40 
41 void RectMatrixCol::Reset (const Matrix& M, int skip, int col, int length)
42 {
43  REPORT
45  ( M.Store()+col+skip*M.Ncols(), length, M.Ncols(), 1 );
46 }
47 
48 void RectMatrixCol::Reset (const Matrix& M, int col)
49 {
50  REPORT
51  RectMatrixRowCol::Reset( M.Store()+col, M.Nrows(), M.Ncols(), 1 );
52 }
53 
54 
56 {
57  REPORT
58  long_Real sum = 0.0; int i = n; Real* s = store; int d = spacing;
59  // while (i--) { sum += (long_Real)*s * *s; s += d; }
60  if (i) for(;;)
61  { sum += (long_Real)*s * *s; if (!(--i)) break; s += d; }
62  return (Real)sum;
63 }
64 
66 {
67  REPORT
68  long_Real sum = 0.0; int i = n;
69  Real* s = store; int d = spacing;
70  Real* s1 = rmrc.store; int d1 = rmrc.spacing;
71  if (i!=rmrc.n)
72  {
73  Tracer tr("newmatrm");
74  Throw(InternalException("Dimensions differ in *"));
75  }
76  // while (i--) { sum += (long_Real)*s * *s1; s += d; s1 += d1; }
77  if (i) for(;;)
78  { sum += (long_Real)*s * *s1; if (!(--i)) break; s += d; s1 += d1; }
79  return (Real)sum;
80 }
81 
83 {
84  REPORT
85  int i = n; Real* s = store; int d = spacing;
86  Real* s1 = rmrc.store; int d1 = rmrc.spacing;
87  if (i!=rmrc.n)
88  {
89  Tracer tr("newmatrm");
90  Throw(InternalException("Dimensions differ in AddScaled"));
91  }
92  // while (i--) { *s += *s1 * r; s += d; s1 += d1; }
93  if (i) for (;;)
94  { *s += *s1 * r; if (!(--i)) break; s += d; s1 += d1; }
95 }
96 
98 {
99  REPORT
100  int i = n; Real* s = store; int d = spacing;
101  Real* s1 = rmrc.store; int d1 = rmrc.spacing;
102  if (i!=rmrc.n)
103  {
104  Tracer tr("newmatrm");
105  Throw(InternalException("Dimensions differ in Divide"));
106  }
107  // while (i--) { *s = *s1 / r; s += d; s1 += d1; }
108  if (i) for (;;) { *s = *s1 / r; if (!(--i)) break; s += d; s1 += d1; }
109 }
110 
112 {
113  REPORT
114  int i = n; Real* s = store; int d = spacing;
115  // while (i--) { *s /= r; s += d; }
116  if (i) for (;;) { *s /= r; if (!(--i)) break; s += d; }
117 }
118 
120 {
121  REPORT
122  int i = n; Real* s = store; int d = spacing;
123  // while (i--) { *s = - *s; s += d; }
124  if (i) for (;;) { *s = - *s; if (!(--i)) break; s += d; }
125 }
126 
128 {
129  REPORT
130  int i = n; Real* s = store; int d = spacing;
131  // while (i--) { *s = 0.0; s += d; }
132  if (i) for (;;) { *s = 0.0; if (!(--i)) break; s += d; }
133 }
134 
136 {
137  REPORT
138  int n = U.n;
139  if (n != V.n)
140  {
141  Tracer tr("newmatrm");
142  Throw(InternalException("Dimensions differ in ComplexScale"));
143  }
144  Real* u = U.store; Real* v = V.store;
145  int su = U.spacing; int sv = V.spacing;
146  //while (n--)
147  //{
148  // Real z = *u * x - *v * y; *v = *u * y + *v * x; *u = z;
149  // u += su; v += sv;
150  //}
151  if (n) for (;;)
152  {
153  Real z = *u * x - *v * y; *v = *u * y + *v * x; *u = z;
154  if (!(--n)) break;
155  u += su; v += sv;
156  }
157 }
158 
160 {
161  REPORT
162  // (U, V) = (U, V) * (c, s) where tau = s/(1+c), c^2 + s^2 = 1
163  int n = U.n;
164  if (n != V.n)
165  {
166  Tracer tr("newmatrm");
167  Throw(InternalException("Dimensions differ in Rotate"));
168  }
169  Real* u = U.store; Real* v = V.store;
170  int su = U.spacing; int sv = V.spacing;
171  //while (n--)
172  //{
173  // Real zu = *u; Real zv = *v;
174  // *u -= s * (zv + zu * tau); *v += s * (zu - zv * tau);
175  // u += su; v += sv;
176  //}
177  if (n) for(;;)
178  {
179  Real zu = *u; Real zv = *v;
180  *u -= s * (zv + zu * tau); *v += s * (zu - zv * tau);
181  if (!(--n)) break;
182  u += su; v += sv;
183  }
184 }
185 
186 
187 // misc procedures for numerical things
188 
189 Real pythag(Real f, Real g, Real& c, Real& s)
190 // return z=sqrt(f*f+g*g), c=f/z, s=g/z
191 // set c=1,s=0 if z==0
192 // avoid floating point overflow or divide by zero
193 {
194  if (f==0 && g==0) { c=1.0; s=0.0; return 0.0; }
195  Real af = f>=0 ? f : -f;
196  Real ag = g>=0 ? g : -g;
197  if (ag<af)
198  {
199  REPORT
200  Real h = g/f; Real sq = sqrt(1.0+h*h);
201  if (f<0) sq = -sq; // make return value non-negative
202  c = 1.0/sq; s = h/sq; return sq*f;
203  }
204  else
205  {
206  REPORT
207  Real h = f/g; Real sq = sqrt(1.0+h*h);
208  if (g<0) sq = -sq;
209  s = 1.0/sq; c = h/sq; return sq*g;
210  }
211 }
212 
213 
214 
215 
216 #ifdef use_namespace
217 }
218 #endif
219 
220 
Real operator*(const RectMatrixRowCol &) const
Definition: newmatrm.cpp:65
void Reset(Real *st, int nx, int sp, int sh)
Definition: newmatrm.h:33
Internal newmat exception - shouldn&#39;t happen.
Definition: newmat.h:2025
void Rotate(RectMatrixCol &U, RectMatrixCol &V, Real tau, Real s)
Definition: newmatrm.cpp:159
double Real
Definition: include.h:307
Real SumSquare() const
Definition: newmatrm.cpp:55
Real pythag(Real f, Real g, Real &c, Real &s)
Definition: newmatrm.cpp:189
int Nrows() const
Definition: newmat.h:494
void Reset(const Matrix &, int, int, int)
Definition: newmatrm.cpp:28
void ComplexScale(RectMatrixCol &U, RectMatrixCol &V, Real x, Real y)
Definition: newmatrm.cpp:135
FloatVector * d
Real * Store() const
Definition: newmat.h:497
void Divide(const RectMatrixRowCol &, Real)
Definition: newmatrm.cpp:97
#define Throw(E)
Definition: myexcept.h:191
The usual rectangular matrix.
Definition: newmat.h:625
int Ncols() const
Definition: newmat.h:495
Real sum(const BaseMatrix &B)
Definition: newmat.h:2105
Real * store
Definition: newmatrm.h:27
void AddScaled(const RectMatrixRowCol &, Real)
Definition: newmatrm.cpp:82
long double long_Real
Definition: include.h:308
#define REPORT
Definition: newmatrm.cpp:21
void Reset(const Matrix &, int, int, int)
Definition: newmatrm.cpp:41


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16