source: ntrip/trunk/BNS/newmat/newmatrm.cpp@ 9042

Last change on this file since 9042 was 810, checked in by mervart, 17 years ago

* empty log message *

File size: 5.1 KB
RevLine 
[810]1/// \ingroup newmat
2///@{
3
4/// \file newmatrm.cpp
5/// Rectangular matrix operations
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
15namespace 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
28void RectMatrixRow::Reset (const Matrix& M, int row, int skip, int length)
29{
30 REPORT
31 RectMatrixRowCol::Reset
32 ( M.Store()+row*M.Ncols()+skip, length, 1, M.Ncols() );
33}
34
35void 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
41void RectMatrixCol::Reset (const Matrix& M, int skip, int col, int length)
42{
43 REPORT
44 RectMatrixRowCol::Reset
45 ( M.Store()+col+skip*M.Ncols(), length, M.Ncols(), 1 );
46}
47
48void RectMatrixCol::Reset (const Matrix& M, int col)
49{
50 REPORT
51 RectMatrixRowCol::Reset( M.Store()+col, M.Nrows(), M.Ncols(), 1 );
52}
53
54
55Real RectMatrixRowCol::SumSquare() const
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
65Real RectMatrixRowCol::operator*(const RectMatrixRowCol& rmrc) const
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
82void RectMatrixRowCol::AddScaled(const RectMatrixRowCol& rmrc, Real r)
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
97void RectMatrixRowCol::Divide(const RectMatrixRowCol& rmrc, Real r)
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
111void RectMatrixRowCol::Divide(Real r)
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
119void RectMatrixRowCol::Negate()
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
127void RectMatrixRowCol::Zero()
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
135void ComplexScale(RectMatrixCol& U, RectMatrixCol& V, Real x, Real y)
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
159void Rotate(RectMatrixCol& U, RectMatrixCol& V, Real tau, Real s)
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
189Real 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
221///@}
Note: See TracBrowser for help on using the repository browser.