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


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