source: nanovis/trunk/newmat11/bandmat.cpp @ 4794

Last change on this file since 4794 was 2096, checked in by ldelgass, 9 years ago

Normalize line endings, set eol-style to native on *.cpp, *.h files

  • Property svn:eol-style set to native
File size: 19.0 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file bandmat.cpp
5/// Band-matrix member functions.
6
7
8// Copyright (C) 1991,2,3,4,9: R B Davies
9
10#define WANT_MATH                    // include.h will get math fns
11
12//#define WANT_STREAM
13
14#include "include.h"
15
16#include "newmat.h"
17#include "newmatrc.h"
18
19#ifdef use_namespace
20namespace NEWMAT {
21#endif
22
23
24
25#ifdef DO_REPORT
26#define REPORT { static ExeCounter ExeCount(__LINE__,10); ++ExeCount; }
27#else
28#define REPORT {}
29#endif
30
31static inline int my_min(int x, int y) { return x < y ? x : y; }
32static inline int my_max(int x, int y) { return x > y ? x : y; }
33
34
35BandMatrix::BandMatrix(const BaseMatrix& M)
36{
37   REPORT // CheckConversion(M);
38   // MatrixConversionCheck mcc;
39   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::BM);
40   GetMatrix(gmx); CornerClear();
41}
42
43void BandMatrix::SetParameters(const GeneralMatrix* gmx)
44{
45   REPORT
46   MatrixBandWidth bw = gmx->bandwidth();
47   lower_val = bw.lower_val; upper_val = bw.upper_val;
48}
49
50void BandMatrix::resize(int n, int lb, int ub)
51{
52   REPORT
53   Tracer tr("BandMatrix::resize");
54   if (lb<0 || ub<0) Throw(ProgramException("Undefined bandwidth"));
55   lower_val = (lb<=n) ? lb : n-1; upper_val = (ub<=n) ? ub : n-1;
56   GeneralMatrix::resize(n,n,n*(lower_val+1+upper_val)); CornerClear();
57}
58
59// SimpleAddOK shows when we can add etc two matrices by a simple vector add
60// and when we can add one matrix into another
61//
62// *gm must be the same type as *this
63// - return 0 if simple add is OK
64// - return 1 if we can add into *gm only
65// - return 2 if we can add into *this only
66// - return 3 if we can't add either way
67//
68// For SP this will still be valid if we swap 1 and 2
69
70/// \brief can we add two band matrices with simple vector add
71///
72/// For band matrices the bandwidths must agree
73
74short BandMatrix::SimpleAddOK(const GeneralMatrix* gm)
75{
76   const BandMatrix* bm = (const BandMatrix*)gm;
77   if (bm->lower_val == lower_val && bm->upper_val == upper_val)
78      { REPORT return 0; }
79   else if (bm->lower_val >= lower_val && bm->upper_val >= upper_val)
80      { REPORT return 1; }
81   else if (bm->lower_val <= lower_val && bm->upper_val <= upper_val)
82      { REPORT return 2; }
83   else { REPORT return 3; }
84}
85
86/// \brief can we add two symmetric band matrices with simple vector add
87///
88/// Sufficient to check lower bandwidths agree
89
90short SymmetricBandMatrix::SimpleAddOK(const GeneralMatrix* gm)
91{
92   const SymmetricBandMatrix* bm = (const SymmetricBandMatrix*)gm;
93   if (bm->lower_val == lower_val) { REPORT return 0; }
94   else if (bm->lower_val > lower_val) { REPORT return 1; }
95   else { REPORT return 2; }
96}
97
98/// \brief resize UpperBandMatrix
99void UpperBandMatrix::resize(int n, int lb, int ub)
100{
101   REPORT
102   if (lb != 0)
103   {
104      Tracer tr("UpperBandMatrix::resize");
105      Throw(ProgramException("UpperBandMatrix with non-zero lower band" ));
106   }
107   BandMatrix::resize(n, lb, ub);
108}
109
110/// \brief resize LowerBandMatrix
111void LowerBandMatrix::resize(int n, int lb, int ub)
112{
113   REPORT
114   if (ub != 0)
115   {
116      Tracer tr("LowerBandMatrix::resize");
117      Throw(ProgramException("LowerBandMatrix with non-zero upper band" ));
118   }
119   BandMatrix::resize(n, lb, ub);
120}
121
122/// \brief resize BandMatrix
123void BandMatrix::resize(const GeneralMatrix& A)
124{
125   REPORT
126   int n = A.Nrows();
127   if (n != A.Ncols())
128   {
129      Tracer tr("BandMatrix::resize(GM)");
130      Throw(NotSquareException(*this));
131   }
132   MatrixBandWidth mbw = A.bandwidth();
133   resize(n, mbw.Lower(), mbw.Upper());
134}
135
136/*
137bool BandMatrix::SameStorageType(const GeneralMatrix& A) const
138{
139   if (type() != A.type()) { REPORT return false; }
140   REPORT
141   return bandwidth() == A.bandwidth();
142}
143
144void BandMatrix::resizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B)
145{
146   REPORT
147   Tracer tr("BandMatrix::resizeForAdd");
148   MatrixBandWidth A_BW = A.bandwidth(); MatrixBandWidth B_BW = B.bandwidth();
149   if ((A_BW.Lower() < 0) | (A_BW.Upper() < 0) | (B_BW.Lower() < 0)
150      | (A_BW.Upper() < 0))
151         Throw(ProgramException("Can't resize to BandMatrix" ));
152   // already know A and B are square
153   resize(A.Nrows(), my_max(A_BW.Lower(), B_BW.Lower()),
154      my_max(A_BW.Upper(), B_BW.Upper()));
155}
156
157void BandMatrix::resizeForSP(const GeneralMatrix& A, const GeneralMatrix& B)
158{
159   REPORT
160   Tracer tr("BandMatrix::resizeForSP");
161   MatrixBandWidth A_BW = A.bandwidth(); MatrixBandWidth B_BW = B.bandwidth();
162   if ((A_BW.Lower() < 0) | (A_BW.Upper() < 0) | (B_BW.Lower() < 0)
163      | (A_BW.Upper() < 0))
164         Throw(ProgramException("Can't resize to BandMatrix" ));
165   // already know A and B are square
166   resize(A.Nrows(), my_min(A_BW.Lower(), B_BW.Lower()),
167      my_min(A_BW.Upper(), B_BW.Upper()));
168}
169*/
170
171/// \brief assignment operator for BandMatrix
172void BandMatrix::operator=(const BaseMatrix& X)
173{
174   REPORT // CheckConversion(X);
175   // MatrixConversionCheck mcc;
176   Eq(X,MatrixType::BM); CornerClear();
177}
178
179/// \brief set unused parts of BandMatrix to zero
180void BandMatrix::CornerClear() const
181{
182   REPORT
183   int i = lower_val; Real* s = store; int bw = lower_val + 1 + upper_val;
184   while (i)
185      { int j = i--; Real* sj = s; s += bw; while (j--) *sj++ = 0.0; }
186   i = upper_val; s = store + storage;
187   while (i)
188      { int j = i--; Real* sj = s; s -= bw; while (j--) *(--sj) = 0.0; }
189}
190
191MatrixBandWidth MatrixBandWidth::operator+(const MatrixBandWidth& bw) const
192{
193   REPORT
194   int l = bw.lower_val; int u = bw.upper_val;
195   l = (lower_val < 0 || l < 0) ? -1 : (lower_val > l) ? lower_val : l;
196   u = (upper_val < 0 || u < 0) ? -1 : (upper_val > u) ? upper_val : u;
197   return MatrixBandWidth(l,u);
198}
199
200MatrixBandWidth MatrixBandWidth::operator*(const MatrixBandWidth& bw) const
201{
202   REPORT
203   int l = bw.lower_val; int u = bw.upper_val;
204   l = (lower_val < 0 || l < 0) ? -1 : lower_val+l;
205   u = (upper_val < 0 || u < 0) ? -1 : upper_val+u;
206   return MatrixBandWidth(l,u);
207}
208
209MatrixBandWidth MatrixBandWidth::minimum(const MatrixBandWidth& bw) const
210{
211   REPORT
212   int l = bw.lower_val; int u = bw.upper_val;
213   if ((lower_val >= 0) && ( (l < 0) || (l > lower_val) )) l = lower_val;
214   if ((upper_val >= 0) && ( (u < 0) || (u > upper_val) )) u = upper_val;
215   return MatrixBandWidth(l,u);
216}
217
218UpperBandMatrix::UpperBandMatrix(const BaseMatrix& M)
219{
220   REPORT // CheckConversion(M);
221   // MatrixConversionCheck mcc;
222   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UB);
223   GetMatrix(gmx); CornerClear();
224}
225
226void UpperBandMatrix::operator=(const BaseMatrix& X)
227{
228   REPORT // CheckConversion(X);
229   // MatrixConversionCheck mcc;
230   Eq(X,MatrixType::UB); CornerClear();
231}
232
233LowerBandMatrix::LowerBandMatrix(const BaseMatrix& M)
234{
235   REPORT // CheckConversion(M);
236   // MatrixConversionCheck mcc;
237   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LB);
238   GetMatrix(gmx); CornerClear();
239}
240
241void LowerBandMatrix::operator=(const BaseMatrix& X)
242{
243   REPORT // CheckConversion(X);
244   // MatrixConversionCheck mcc;
245   Eq(X,MatrixType::LB); CornerClear();
246}
247
248BandLUMatrix::BandLUMatrix(const BaseMatrix& m)
249{
250   REPORT
251   Tracer tr("BandLUMatrix");
252   storage2 = 0; store2 = 0; indx = 0; // in event of exception during build
253   GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate();
254   if (gm->nrows() != gm->ncols())
255      { gm->tDelete(); Throw(NotSquareException(*this)); }
256   if (gm->type() == MatrixType::BC)
257      { REPORT  ((BandLUMatrix*)gm)->get_aux(*this); GetMatrix(gm); }
258   else
259   {
260      REPORT
261      BandMatrix* gm1 = (BandMatrix*)(gm->Evaluate(MatrixType::BM));
262      m1 = gm1->lower_val; m2 = gm1->upper_val;
263      GetMatrix(gm1);
264      d = true; sing = false;
265      indx = new int [nrows_val]; MatrixErrorNoSpace(indx);
266      MONITOR_INT_NEW("Index (BndLUMat)",nrows_val,indx)
267      storage2 = nrows_val * m1;
268      store2 = new Real [storage2]; MatrixErrorNoSpace(store2);
269      MONITOR_REAL_NEW("Make (BandLUMat)",storage2,store2)
270      ludcmp();
271   }
272}
273
274GeneralMatrix* BandLUMatrix::Evaluate(MatrixType mt)
275{
276   if (Compare(this->Type(),mt)) { REPORT return this; }
277   REPORT
278   Tracer et("BandLUMatrix::Evaluate");
279   bool dummy = true;
280   if (dummy) Throw(ProgramException("Illegal use of BandLUMatrix", *this));
281   return this;
282}
283
284// could we use SetParameters instead of this
285void BandLUMatrix::get_aux(BandLUMatrix& X)
286{
287   X.d = d; X.sing = sing; X.storage2 = storage2; X.m1 = m1; X.m2 = m2;   
288   if (tag_val == 0 || tag_val == 1) // reuse the array
289   {
290      REPORT
291      X.indx = indx; indx = 0;
292      X.store2 = store2; store2 = 0;
293      d = true; sing = true; storage2 = 0; m1 = 0; m2 = 0;
294      return;
295   }
296   else if (nrows_val == 0)
297   {
298      REPORT
299      indx = 0; store2 = 0; storage2 = 0;
300      d = true; sing = true; m1 = m2 = 0;
301      return;
302   }
303   else                              // copy the array
304   {
305      REPORT
306      Tracer tr("BandLUMatrix::get_aux");
307      int *ix = new int [nrows_val]; MatrixErrorNoSpace(ix);
308      MONITOR_INT_NEW("Index (BLUM::get_aux)", nrows_val, ix)
309      int n = nrows_val; int* i = ix; int* j = indx;
310      while(n--) *i++ = *j++;
311      X.indx = ix;
312      Real *rx = new Real [storage2]; MatrixErrorNoSpace(indx);
313      MONITOR_REAL_NEW("Index (BLUM::get_aux)", storage2, rx)
314      newmat_block_copy(storage2, store2, rx);
315      X.store2 = rx;
316   }
317}
318
319BandLUMatrix::BandLUMatrix(const BandLUMatrix& gm) : GeneralMatrix()
320{
321   REPORT
322   Tracer tr("BandLUMatrix(const BandLUMatrix&)");
323   ((BandLUMatrix&)gm).get_aux(*this);
324   GetMatrix(&gm);
325}
326
327void BandLUMatrix::operator=(const BandLUMatrix& gm)
328{
329   if (&gm == this) { REPORT tag_val = -1; return; }
330   REPORT
331   delete [] indx; indx = 0;
332   delete [] store2; store2 = 0; storage2 = 0;
333   ((BandLUMatrix&)gm).get_aux(*this);
334   Eq(gm);
335}
336   
337
338
339
340
341
342
343
344BandLUMatrix::~BandLUMatrix()
345{
346   REPORT
347   MONITOR_INT_DELETE("Index (BndLUMat)",nrows_val,indx)
348   MONITOR_REAL_DELETE("Delete (BndLUMt)",storage2,store2)
349   delete [] indx; delete [] store2;
350}
351
352MatrixType BandLUMatrix::type() const { REPORT return MatrixType::BC; }
353
354
355LogAndSign BandLUMatrix::log_determinant() const
356{
357   REPORT
358   if (sing) return 0.0;
359   Real* a = store; int w = m1+1+m2; LogAndSign sum; int i = nrows_val;
360   // while (i--) { sum *= *a; a += w; }
361   if (i) for (;;) { sum *= *a; if (!(--i)) break; a += w; }
362   if (!d) sum.ChangeSign(); return sum;
363}
364
365GeneralMatrix* BandMatrix::MakeSolver()
366{
367   REPORT
368   GeneralMatrix* gm = new BandLUMatrix(*this);
369   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
370}
371
372
373void BandLUMatrix::ludcmp()
374{
375   REPORT
376   Real* a = store2; int i = storage2;
377   // clear store2 - so unused locations are always zero -
378   // required by operator==
379   while (i--) *a++ = 0.0;
380   a = store;
381   i = m1; int j = m2; int k; int n = nrows_val; int w = m1 + 1 + m2;
382   while (i)
383   {
384      Real* ai = a + i;
385      k = ++j; while (k--) *a++ = *ai++;
386      k = i--; while (k--) *a++ = 0.0;
387   }
388
389   a = store; int l = m1;
390   for (k=0; k<n; k++)
391   {
392      Real x = *a; i = k; Real* aj = a;
393      if (l < n) l++;
394      for (j=k+1; j<l; j++)
395         { aj += w; if (fabs(x) < fabs(*aj)) { x = *aj; i = j; } }
396      indx[k] = i;
397      if (x==0) { sing = true; return; }
398      if (i!=k)
399      {
400         d = !d; Real* ak = a; Real* ai = store + i * w; j = w;
401         while (j--) { x = *ak; *ak++ = *ai; *ai++ = x; }
402      }
403      aj = a + w; Real* m = store2 + m1 * k;
404      for (j=k+1; j<l; j++)
405      {
406         *m++ = x = *aj / *a; i = w; Real* ak = a;
407         while (--i) { Real* aj1 = aj++; *aj1 = *aj - x * *(++ak); }
408         *aj++ = 0.0;
409      }
410      a += w;
411   }
412}
413
414void BandLUMatrix::lubksb(Real* B, int mini)
415{
416   REPORT
417   Tracer tr("BandLUMatrix::lubksb");
418   if (sing) Throw(SingularException(*this));
419   int n = nrows_val; int l = m1; int w = m1 + 1 + m2;
420
421   for (int k=0; k<n; k++)
422   {
423      int i = indx[k];
424      if (i!=k) { Real x=B[k]; B[k]=B[i]; B[i]=x; }
425      if (l<n) l++;
426      Real* m = store2 + k*m1; Real* b = B+k; Real* bi = b;
427      for (i=k+1; i<l; i++)  *(++bi) -= *m++ * *b;
428   }
429
430   l = -m1;
431   for (int i = n-1; i>=mini; i--)
432   {
433      Real* b = B + i; Real* bk = b; Real x = *bk;
434      Real* a = store + w*i; Real y = *a;
435      int k = l+m1; while (k--) x -=  *(++a) * *(++bk);
436      *b = x / y;
437      if (l < m2) l++;
438   }
439}
440
441void BandLUMatrix::Solver(MatrixColX& mcout, const MatrixColX& mcin)
442{
443   REPORT
444   int i = mcin.skip; Real* el = mcin.data-i; Real* el1=el;
445   while (i--) *el++ = 0.0;
446   el += mcin.storage; i = nrows_val - mcin.skip - mcin.storage;
447   while (i--) *el++ = 0.0;
448   lubksb(el1, mcout.skip);
449}
450
451// Do we need check for entirely zero output?
452
453
454void UpperBandMatrix::Solver(MatrixColX& mcout,
455   const MatrixColX& mcin)
456{
457   REPORT
458   int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
459   while (i-- > 0) *elx++ = 0.0;
460   int nr = mcin.skip+mcin.storage;
461   elx = mcin.data+mcin.storage; Real* el = elx;
462   int j = mcout.skip+mcout.storage-nr; i = nr-mcout.skip;
463   while (j-- > 0) *elx++ = 0.0;
464
465   Real* Ael = store + (upper_val+1)*(i-1)+1; j = 0;
466   if (i > 0) for(;;)
467   {
468      elx = el; Real sum = 0.0; int jx = j;
469      while (jx--) sum += *(--Ael) * *(--elx);
470      elx--; *elx = (*elx - sum) / *(--Ael);
471      if (--i <= 0) break;
472      if (j<upper_val) Ael -= upper_val - (++j); else el--;
473   }
474}
475
476void LowerBandMatrix::Solver(MatrixColX& mcout,
477   const MatrixColX& mcin)
478{
479   REPORT
480   int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
481   while (i-- > 0) *elx++ = 0.0;
482   int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.data+mcin.storage;
483   int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
484   while (j-- > 0) *elx++ = 0.0;
485
486   Real* el = mcin.data;
487   Real* Ael = store + (lower_val+1)*nc + lower_val;
488   j = 0;
489   if (i > 0) for(;;)
490   {
491      elx = el; Real sum = 0.0; int jx = j;
492      while (jx--) sum += *Ael++ * *elx++;
493      *elx = (*elx - sum) / *Ael++;
494      if (--i <= 0) break;
495      if (j<lower_val) Ael += lower_val - (++j); else el++;
496   }
497}
498
499
500LogAndSign BandMatrix::log_determinant() const
501{
502   REPORT
503   BandLUMatrix C(*this); return C.log_determinant();
504}
505
506LogAndSign LowerBandMatrix::log_determinant() const
507{
508   REPORT
509   int i = nrows_val; LogAndSign sum;
510   Real* s = store + lower_val; int j = lower_val + 1;
511//   while (i--) { sum *= *s; s += j; }
512   if (i) for (;;) { sum *= *s; if (!(--i)) break; s += j; }
513   ((GeneralMatrix&)*this).tDelete(); return sum;
514}
515
516LogAndSign UpperBandMatrix::log_determinant() const
517{
518   REPORT
519   int i = nrows_val; LogAndSign sum; Real* s = store; int j = upper_val + 1;
520//   while (i--) { sum *= *s; s += j; }
521   if (i) for (;;) { sum *= *s; if (!(--i)) break; s += j; }
522   ((GeneralMatrix&)*this).tDelete(); return sum;
523}
524
525GeneralMatrix* SymmetricBandMatrix::MakeSolver()
526{
527   REPORT
528   GeneralMatrix* gm = new BandLUMatrix(*this);
529   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
530}
531
532SymmetricBandMatrix::SymmetricBandMatrix(const BaseMatrix& M)
533{
534   REPORT  // CheckConversion(M);
535   // MatrixConversionCheck mcc;
536   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::SB);
537   GetMatrix(gmx);
538}
539
540GeneralMatrix* SymmetricBandMatrix::Transpose(TransposedMatrix*, MatrixType mt)
541{ REPORT  return Evaluate(mt); }
542
543LogAndSign SymmetricBandMatrix::log_determinant() const
544{
545   REPORT
546   BandLUMatrix C(*this); return C.log_determinant();
547}
548
549void SymmetricBandMatrix::SetParameters(const GeneralMatrix* gmx)
550{ REPORT lower_val = gmx->bandwidth().lower_val; }
551
552void SymmetricBandMatrix::resize(int n, int lb)
553{
554   REPORT
555   Tracer tr("SymmetricBandMatrix::resize");
556   if (lb<0) Throw(ProgramException("Undefined bandwidth"));
557   lower_val = (lb<=n) ? lb : n-1;
558   GeneralMatrix::resize(n,n,n*(lower_val+1));
559}
560
561void SymmetricBandMatrix::resize(const GeneralMatrix& A)
562{
563   REPORT
564   int n = A.Nrows();
565   if (n != A.Ncols())
566   {
567      Tracer tr("SymmetricBandMatrix::resize(GM)");
568      Throw(NotSquareException(*this));
569   }
570   MatrixBandWidth mbw = A.bandwidth(); int b = mbw.Lower();
571   if (b != mbw.Upper())
572   {
573      Tracer tr("SymmetricBandMatrix::resize(GM)");
574      Throw(ProgramException("Upper and lower band-widths not equal"));
575   }
576   resize(n, b);
577}
578/*
579bool SymmetricBandMatrix::SameStorageType(const GeneralMatrix& A) const
580{
581   if (type() != A.type()) { REPORT return false; }
582   REPORT
583   return bandwidth() == A.bandwidth();
584}
585
586void SymmetricBandMatrix::resizeForAdd(const GeneralMatrix& A,
587   const GeneralMatrix& B)
588{
589   REPORT
590   Tracer tr("SymmetricBandMatrix::resizeForAdd");
591   MatrixBandWidth A_BW = A.bandwidth(); MatrixBandWidth B_BW = B.bandwidth();
592   if ((A_BW.Lower() < 0) | (B_BW.Lower() < 0))
593         Throw(ProgramException("Can't resize to SymmetricBandMatrix" ));
594   // already know A and B are square
595   resize(A.Nrows(), my_max(A_BW.Lower(), B_BW.Lower()));
596}
597
598void SymmetricBandMatrix::resizeForSP(const GeneralMatrix& A,
599   const GeneralMatrix& B)
600{
601   REPORT
602   Tracer tr("SymmetricBandMatrix::resizeForSP");
603   MatrixBandWidth A_BW = A.bandwidth(); MatrixBandWidth B_BW = B.bandwidth();
604   if ((A_BW.Lower() < 0) | (B_BW.Lower() < 0))
605         Throw(ProgramException("Can't resize to SymmetricBandMatrix" ));
606   // already know A and B are square
607   resize(A.Nrows(), my_min(A_BW.Lower(), B_BW.Lower()));
608}
609*/
610
611void SymmetricBandMatrix::operator=(const BaseMatrix& X)
612{
613   REPORT // CheckConversion(X);
614   // MatrixConversionCheck mcc;
615   Eq(X,MatrixType::SB);
616}
617
618void SymmetricBandMatrix::CornerClear() const
619{
620   // set unused parts of BandMatrix to zero
621   REPORT
622   int i = lower_val; Real* s = store; int bw = lower_val + 1;
623   if (i) for(;;)
624   {
625      int j = i;
626      Real* sj = s;
627      while (j--) *sj++ = 0.0;
628      if (!(--i)) break;
629      s += bw;
630   }
631}
632
633MatrixBandWidth SymmetricBandMatrix::bandwidth() const
634   { REPORT return MatrixBandWidth(lower_val,lower_val); }
635
636GeneralMatrix* BandMatrix::Image() const
637{
638   REPORT
639   GeneralMatrix* gm = new BandMatrix(*this); MatrixErrorNoSpace(gm);
640   return gm;
641}
642
643GeneralMatrix* UpperBandMatrix::Image() const
644{
645   REPORT
646   GeneralMatrix* gm = new UpperBandMatrix(*this); MatrixErrorNoSpace(gm);
647   return gm;
648}
649
650GeneralMatrix* LowerBandMatrix::Image() const
651{
652   REPORT
653   GeneralMatrix* gm = new LowerBandMatrix(*this); MatrixErrorNoSpace(gm);
654   return gm;
655}
656
657GeneralMatrix* SymmetricBandMatrix::Image() const
658{
659   REPORT
660   GeneralMatrix* gm = new SymmetricBandMatrix(*this); MatrixErrorNoSpace(gm);
661   return gm;
662}
663
664GeneralMatrix* BandLUMatrix::Image() const
665{
666   REPORT
667   GeneralMatrix* gm = new BandLUMatrix(*this); MatrixErrorNoSpace(gm);
668   return gm;
669}
670
671
672inline Real square(Real x) { return x*x; }
673
674Real SymmetricBandMatrix::sum_square() const
675{
676   REPORT
677   CornerClear();
678   Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows_val;
679   int l=lower_val;
680   while (i--)
681      { int j = l; while (j--) sum2 += square(*s++); sum1 += square(*s++); }
682   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
683}
684
685Real SymmetricBandMatrix::sum_absolute_value() const
686{
687   REPORT
688   CornerClear();
689   Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows_val;
690   int l=lower_val;
691   while (i--)
692      { int j = l; while (j--) sum2 += fabs(*s++); sum1 += fabs(*s++); }
693   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
694}
695
696Real SymmetricBandMatrix::sum() const
697{
698   REPORT
699   CornerClear();
700   Real sum1=0.0; Real sum2=0.0; Real* s=store; int i=nrows_val;
701   int l=lower_val;
702   while (i--)
703      { int j = l; while (j--) sum2 += *s++; sum1 += *s++; }
704   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
705}
706
707
708
709
710
711#ifdef use_namespace
712}
713#endif
714
715///@}
716
717
Note: See TracBrowser for help on using the repository browser.