source: ntrip/trunk/BNC/newmat/newmat5.cpp@ 2013

Last change on this file since 2013 was 2013, checked in by mervart, 13 years ago

* empty log message *

File size: 14.7 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file newmat5.cpp
5/// Transpose, evaluate, operations with scalar, matrix input.
6
7// Copyright (C) 1991,2,3,4: R B Davies
8
9//#define WANT_STREAM
10
11#include "include.h"
12
13#include "newmat.h"
14#include "newmatrc.h"
15
16#ifdef use_namespace
17namespace NEWMAT {
18#endif
19
20
21#ifdef DO_REPORT
22#define REPORT { static ExeCounter ExeCount(__LINE__,5); ++ExeCount; }
23#else
24#define REPORT {}
25#endif
26
27
28/************************ carry out operations ******************************/
29
30
31GeneralMatrix* GeneralMatrix::Transpose(TransposedMatrix* tm, MatrixType mt)
32{
33 GeneralMatrix* gm1;
34
35 if (Compare(Type().t(),mt))
36 {
37 REPORT
38 gm1 = mt.New(ncols_val,nrows_val,tm);
39 for (int i=0; i<ncols_val; i++)
40 {
41 MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
42 MatrixCol mc(this, mr.Data(), LoadOnEntry, i);
43 }
44 }
45 else
46 {
47 REPORT
48 gm1 = mt.New(ncols_val,nrows_val,tm);
49 MatrixRow mr(this, LoadOnEntry);
50 MatrixCol mc(gm1, StoreOnExit+DirectPart);
51 int i = nrows_val;
52 while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
53 }
54 tDelete(); gm1->ReleaseAndDelete(); return gm1;
55}
56
57GeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
58{ REPORT return Evaluate(mt); }
59
60
61GeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
62{ REPORT return Evaluate(mt); }
63
64GeneralMatrix* ColumnVector::Transpose(TransposedMatrix*, MatrixType mt)
65{
66 REPORT
67 GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
68 gmx->nrows_val = 1; gmx->ncols_val = gmx->storage = storage;
69 return BorrowStore(gmx,mt);
70}
71
72GeneralMatrix* RowVector::Transpose(TransposedMatrix*, MatrixType mt)
73{
74 REPORT
75 GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
76 gmx->ncols_val = 1; gmx->nrows_val = gmx->storage = storage;
77 return BorrowStore(gmx,mt);
78}
79
80GeneralMatrix* IdentityMatrix::Transpose(TransposedMatrix*, MatrixType mt)
81{ REPORT return Evaluate(mt); }
82
83GeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
84{
85 if (Compare(this->Type(),mt)) { REPORT return this; }
86 REPORT
87 GeneralMatrix* gmx = mt.New(nrows_val,ncols_val,this);
88 MatrixRow mr(this, LoadOnEntry);
89 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
90 int i=nrows_val;
91 while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
92 tDelete(); gmx->ReleaseAndDelete(); return gmx;
93}
94
95GeneralMatrix* CroutMatrix::Evaluate(MatrixType mt)
96{
97 if (Compare(this->Type(),mt)) { REPORT return this; }
98 REPORT
99 Tracer et("CroutMatrix::Evaluate");
100 bool dummy = true;
101 if (dummy) Throw(ProgramException("Illegal use of CroutMatrix", *this));
102 return this;
103}
104
105GeneralMatrix* GenericMatrix::Evaluate(MatrixType mt)
106 { REPORT return gm->Evaluate(mt); }
107
108GeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
109{
110 gm=((BaseMatrix*&)bm)->Evaluate();
111 int nr=gm->Nrows(); int nc=gm->Ncols();
112 Compare(gm->Type().AddEqualEl(),mt);
113 if (!(mt==gm->Type()))
114 {
115 REPORT
116 GeneralMatrix* gmx = mt.New(nr,nc,this);
117 MatrixRow mr(gm, LoadOnEntry);
118 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
119 while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
120 gmx->ReleaseAndDelete(); gm->tDelete();
121 return gmx;
122 }
123 else if (gm->reuse())
124 {
125 REPORT gm->Add(f);
126 return gm;
127 }
128 else
129 {
130 REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
131 gmy->ReleaseAndDelete(); gmy->Add(gm,f);
132 return gmy;
133 }
134}
135
136GeneralMatrix* NegShiftedMatrix::Evaluate(MatrixType mt)
137{
138 gm=((BaseMatrix*&)bm)->Evaluate();
139 int nr=gm->Nrows(); int nc=gm->Ncols();
140 Compare(gm->Type().AddEqualEl(),mt);
141 if (!(mt==gm->Type()))
142 {
143 REPORT
144 GeneralMatrix* gmx = mt.New(nr,nc,this);
145 MatrixRow mr(gm, LoadOnEntry);
146 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
147 while (nr--) { mrx.NegAdd(mr,f); mrx.Next(); mr.Next(); }
148 gmx->ReleaseAndDelete(); gm->tDelete();
149 return gmx;
150 }
151 else if (gm->reuse())
152 {
153 REPORT gm->NegAdd(f);
154 return gm;
155 }
156 else
157 {
158 REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
159 gmy->ReleaseAndDelete(); gmy->NegAdd(gm,f);
160 return gmy;
161 }
162}
163
164GeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
165{
166 gm=((BaseMatrix*&)bm)->Evaluate();
167 int nr=gm->Nrows(); int nc=gm->Ncols();
168 if (Compare(gm->Type(),mt))
169 {
170 if (gm->reuse())
171 {
172 REPORT gm->Multiply(f);
173 return gm;
174 }
175 else
176 {
177 REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
178 gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
179 return gmx;
180 }
181 }
182 else
183 {
184 REPORT
185 GeneralMatrix* gmx = mt.New(nr,nc,this);
186 MatrixRow mr(gm, LoadOnEntry);
187 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
188 while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
189 gmx->ReleaseAndDelete(); gm->tDelete();
190 return gmx;
191 }
192}
193
194GeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
195{
196 gm=((BaseMatrix*&)bm)->Evaluate();
197 int nr=gm->Nrows(); int nc=gm->Ncols();
198 if (Compare(gm->Type(),mt))
199 {
200 if (gm->reuse())
201 {
202 REPORT gm->Negate();
203 return gm;
204 }
205 else
206 {
207 REPORT
208 GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
209 gmx->ReleaseAndDelete(); gmx->Negate(gm);
210 return gmx;
211 }
212 }
213 else
214 {
215 REPORT
216 GeneralMatrix* gmx = mt.New(nr,nc,this);
217 MatrixRow mr(gm, LoadOnEntry);
218 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
219 while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
220 gmx->ReleaseAndDelete(); gm->tDelete();
221 return gmx;
222 }
223}
224
225GeneralMatrix* ReversedMatrix::Evaluate(MatrixType mt)
226{
227 gm=((BaseMatrix*&)bm)->Evaluate(); GeneralMatrix* gmx;
228
229 if ((gm->Type()).is_band() && ! (gm->Type()).is_diagonal())
230 {
231 gm->tDelete();
232 Throw(NotDefinedException("Reverse", "band matrices"));
233 }
234
235 if (gm->reuse()) { REPORT gm->ReverseElements(); gmx = gm; }
236 else
237 {
238 REPORT
239 gmx = gm->Type().New(gm->Nrows(), gm->Ncols(), this);
240 gmx->ReverseElements(gm); gmx->ReleaseAndDelete();
241 }
242 return gmx->Evaluate(mt); // target matrix is different type?
243
244}
245
246GeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
247{
248 REPORT
249 gm=((BaseMatrix*&)bm)->Evaluate();
250 Compare(gm->Type().t(),mt);
251 GeneralMatrix* gmx=gm->Transpose(this, mt);
252 return gmx;
253}
254
255GeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
256{
257 gm = ((BaseMatrix*&)bm)->Evaluate();
258 GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
259 gmx->nrows_val = 1; gmx->ncols_val = gmx->storage = gm->storage;
260 return gm->BorrowStore(gmx,mt);
261}
262
263GeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
264{
265 gm = ((BaseMatrix*&)bm)->Evaluate();
266 GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
267 gmx->ncols_val = 1; gmx->nrows_val = gmx->storage = gm->storage;
268 return gm->BorrowStore(gmx,mt);
269}
270
271GeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
272{
273 gm = ((BaseMatrix*&)bm)->Evaluate();
274 GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
275 gmx->nrows_val = gmx->ncols_val = gmx->storage = gm->storage;
276 return gm->BorrowStore(gmx,mt);
277}
278
279GeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
280{
281 Tracer tr("MatedMatrix::Evaluate");
282 gm = ((BaseMatrix*&)bm)->Evaluate();
283 GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
284 gmx->nrows_val = nr; gmx->ncols_val = nc; gmx->storage = gm->storage;
285 if (nr*nc != gmx->storage)
286 Throw(IncompatibleDimensionsException());
287 return gm->BorrowStore(gmx,mt);
288}
289
290GeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
291{
292 REPORT
293 Tracer tr("SubMatrix(evaluate)");
294 gm = ((BaseMatrix*&)bm)->Evaluate();
295 if (row_number < 0) row_number = gm->Nrows();
296 if (col_number < 0) col_number = gm->Ncols();
297 if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
298 {
299 gm->tDelete();
300 Throw(SubMatrixDimensionException());
301 }
302 if (IsSym) Compare(gm->Type().ssub(), mt);
303 else Compare(gm->Type().sub(), mt);
304 GeneralMatrix* gmx = mt.New(row_number, col_number, this);
305 int i = row_number;
306 MatrixRow mr(gm, LoadOnEntry, row_skip);
307 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
308 MatrixRowCol sub;
309 while (i--)
310 {
311 mr.SubRowCol(sub, col_skip, col_number); // put values in sub
312 mrx.Copy(sub); mrx.Next(); mr.Next();
313 }
314 gmx->ReleaseAndDelete(); gm->tDelete();
315 return gmx;
316}
317
318
319GeneralMatrix* ReturnMatrix::Evaluate(MatrixType mt)
320{
321 return gm->Evaluate(mt);
322}
323
324
325void GeneralMatrix::Add(GeneralMatrix* gm1, Real f)
326{
327 REPORT
328 Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
329 while (i--)
330 { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; }
331 i = storage & 3; while (i--) *s++ = *s1++ + f;
332}
333
334void GeneralMatrix::Add(Real f)
335{
336 REPORT
337 Real* s=store; int i=(storage >> 2);
338 while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; }
339 i = storage & 3; while (i--) *s++ += f;
340}
341
342void GeneralMatrix::NegAdd(GeneralMatrix* gm1, Real f)
343{
344 REPORT
345 Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
346 while (i--)
347 { *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; }
348 i = storage & 3; while (i--) *s++ = f - *s1++;
349}
350
351void GeneralMatrix::NegAdd(Real f)
352{
353 REPORT
354 Real* s=store; int i=(storage >> 2);
355 while (i--)
356 {
357 *s = f - *s; s++; *s = f - *s; s++;
358 *s = f - *s; s++; *s = f - *s; s++;
359 }
360 i = storage & 3; while (i--) { *s = f - *s; s++; }
361}
362
363void GeneralMatrix::Negate(GeneralMatrix* gm1)
364{
365 // change sign of elements
366 REPORT
367 Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
368 while (i--)
369 { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); }
370 i = storage & 3; while(i--) *s++ = -(*s1++);
371}
372
373void GeneralMatrix::Negate()
374{
375 REPORT
376 Real* s=store; int i=(storage >> 2);
377 while (i--)
378 { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; }
379 i = storage & 3; while(i--) { *s = -(*s); s++; }
380}
381
382void GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f)
383{
384 REPORT
385 Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
386 while (i--)
387 { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; }
388 i = storage & 3; while (i--) *s++ = *s1++ * f;
389}
390
391void GeneralMatrix::Multiply(Real f)
392{
393 REPORT
394 Real* s=store; int i=(storage >> 2);
395 while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; }
396 i = storage & 3; while (i--) *s++ *= f;
397}
398
399
400/************************ MatrixInput routines ****************************/
401
402// int MatrixInput::n; // number values still to be read
403// Real* MatrixInput::r; // pointer to next location to be read to
404
405MatrixInput MatrixInput::operator<<(double f)
406{
407 REPORT
408 Tracer et("MatrixInput");
409 if (n<=0) Throw(ProgramException("List of values too long"));
410 *r = (Real)f; int n1 = n-1; n=0; // n=0 so we won't trigger exception
411 return MatrixInput(n1, r+1);
412}
413
414
415MatrixInput GeneralMatrix::operator<<(double f)
416{
417 REPORT
418 Tracer et("MatrixInput");
419 int n = Storage();
420 if (n<=0) Throw(ProgramException("Loading data to zero length matrix"));
421 Real* r; r = Store(); *r = (Real)f; n--;
422 return MatrixInput(n, r+1);
423}
424
425MatrixInput GetSubMatrix::operator<<(double f)
426{
427 REPORT
428 Tracer et("MatrixInput (GetSubMatrix)");
429 SetUpLHS();
430 if (row_number != 1 || col_skip != 0 || col_number != gm->Ncols())
431 {
432 Throw(ProgramException("MatrixInput requires complete rows"));
433 }
434 MatrixRow mr(gm, DirectPart, row_skip); // to pick up location and length
435 int n = mr.Storage();
436 if (n<=0)
437 {
438 Throw(ProgramException("Loading data to zero length row"));
439 }
440 Real* r; r = mr.Data(); *r = (Real)f; n--;
441 if (+(mr.cw*HaveStore))
442 {
443 Throw(ProgramException("Fails with this matrix type"));
444 }
445 return MatrixInput(n, r+1);
446}
447
448MatrixInput MatrixInput::operator<<(float f)
449{
450 REPORT
451 Tracer et("MatrixInput");
452 if (n<=0) Throw(ProgramException("List of values too long"));
453 *r = (Real)f; int n1 = n-1; n=0; // n=0 so we won't trigger exception
454 return MatrixInput(n1, r+1);
455}
456
457
458MatrixInput GeneralMatrix::operator<<(float f)
459{
460 REPORT
461 Tracer et("MatrixInput");
462 int n = Storage();
463 if (n<=0) Throw(ProgramException("Loading data to zero length matrix"));
464 Real* r; r = Store(); *r = (Real)f; n--;
465 return MatrixInput(n, r+1);
466}
467
468MatrixInput GetSubMatrix::operator<<(float f)
469{
470 REPORT
471 Tracer et("MatrixInput (GetSubMatrix)");
472 SetUpLHS();
473 if (row_number != 1 || col_skip != 0 || col_number != gm->Ncols())
474 {
475 Throw(ProgramException("MatrixInput requires complete rows"));
476 }
477 MatrixRow mr(gm, DirectPart, row_skip); // to pick up location and length
478 int n = mr.Storage();
479 if (n<=0)
480 {
481 Throw(ProgramException("Loading data to zero length row"));
482 }
483 Real* r; r = mr.Data(); *r = (Real)f; n--;
484 if (+(mr.cw*HaveStore))
485 {
486 Throw(ProgramException("Fails with this matrix type"));
487 }
488 return MatrixInput(n, r+1);
489}
490MatrixInput::~MatrixInput()
491{
492 REPORT
493 Tracer et("MatrixInput");
494 if (n!=0) Throw(ProgramException("A list of values was too short"));
495}
496
497MatrixInput BandMatrix::operator<<(double)
498{
499 Tracer et("MatrixInput");
500 bool dummy = true;
501 if (dummy) // get rid of warning message
502 Throw(ProgramException("Cannot use list read with a BandMatrix"));
503 return MatrixInput(0, 0);
504}
505
506MatrixInput BandMatrix::operator<<(float)
507{
508 Tracer et("MatrixInput");
509 bool dummy = true;
510 if (dummy) // get rid of warning message
511 Throw(ProgramException("Cannot use list read with a BandMatrix"));
512 return MatrixInput(0, 0);
513}
514
515void BandMatrix::operator<<(const double*)
516{ Throw(ProgramException("Cannot use array read with a BandMatrix")); }
517
518void BandMatrix::operator<<(const float*)
519{ Throw(ProgramException("Cannot use array read with a BandMatrix")); }
520
521void BandMatrix::operator<<(const int*)
522{ Throw(ProgramException("Cannot use array read with a BandMatrix")); }
523
524void SymmetricBandMatrix::operator<<(const double*)
525{ Throw(ProgramException("Cannot use array read with a BandMatrix")); }
526
527void SymmetricBandMatrix::operator<<(const float*)
528{ Throw(ProgramException("Cannot use array read with a BandMatrix")); }
529
530void SymmetricBandMatrix::operator<<(const int*)
531{ Throw(ProgramException("Cannot use array read with a BandMatrix")); }
532
533// ************************* Reverse order of elements ***********************
534
535void GeneralMatrix::ReverseElements(GeneralMatrix* gm)
536{
537 // reversing into a new matrix
538 REPORT
539 int n = Storage(); Real* rx = Store() + n; Real* x = gm->Store();
540 while (n--) *(--rx) = *(x++);
541}
542
543void GeneralMatrix::ReverseElements()
544{
545 // reversing in place
546 REPORT
547 int n = Storage(); Real* x = Store(); Real* rx = x + n;
548 n /= 2;
549 while (n--) { Real t = *(--rx); *rx = *x; *(x++) = t; }
550}
551
552
553#ifdef use_namespace
554}
555#endif
556
557///@}
Note: See TracBrowser for help on using the repository browser.