Toast++  1.0.2 (r.539)
Forward and inverse modelling in optical tomography
matrix.h
1 // -*-C++-*-
2 // ==========================================================================
3 // Module mathlib
4 // File matrix.h
5 // Declaration of template class TMatrix
6 // This is the root for all matrix classes and defines generic matrix methods
7 //
8 // Inheritance:
9 // ------------
10 // TMatrix ----> ...
11 // ==========================================================================
12 
13 #ifndef __MATRIX_H
14 #define __MATRIX_H
15 
16 #include "vector.h"
17 
19 
20 enum MatrixStorage {
21  MATRIX_DENSE,
22  MATRIX_SYM,
23  MATRIX_DIAG,
24  MATRIX_COORD,
25  MATRIX_COMPROW,
26  MATRIX_COMPCOL,
27  MATRIX_SYMCOMPROW
28 };
30 
31 template<class MT> class TSymMatrix;
32 template<class MT> class TPreconditioner;
33 
34 struct IterativeSolverResult {
35  int it_count;
36  double rel_err;
37 };
38 
39 // ==========================================================================
40 // Nonmember declarations
41 // ==========================================================================
42 
43 template<class MT> class TMatrix;
44 
45 template<class MT>
46 TVector<MT> Ax (const TMatrix<MT> &A, const TVector<MT> &x)
47 { TVector<MT> b; A.Ax(x,b); return b; }
48 // Returns the product A x
49 
50 template<class MT>
51 TVector<MT> ATx (const TMatrix<MT> &A, const TVector<MT> &x)
52 { TVector<MT> b; A.ATx(x,b); return b; }
53 // Returns the product A^T x
54 
55 template<class MT>
56 TSymMatrix<MT> ATA (const TMatrix<MT> &A);
57 // Returns symmetric matrix A^T A (dimension A.Cols() x A.Cols())
58 
59 template<class MT>
60 TSymMatrix<MT> AAT (const TMatrix<MT> &A);
61 // Returns symmetric matrix A A^T (dimension A.Rows() x A.Rows())
62 
63 template<class MT>
64 TVector<MT> ATA_diag (const TMatrix<MT> &A);
65 // Returns the diagonal of A^T A as a vector
66 
67 template<class MT>
68 int PCG (const TMatrix<MT> &A, const TVector<MT> &b, TVector<MT> &x,
69  double &tol, TPreconditioner<MT> *precon = 0, int maxit = 0);
70 // PCG (preconditioned conjugate gradient) linear solver
71 // Solves Ax = b for given tolerance and preconditioner
72 
73 template<class MT>
74 void PCG (const TMatrix<MT> &A, const TVector<MT> *b, TVector<MT> *x,
75  int nrhs, double tol, int maxit = 0, TPreconditioner<MT> *precon = 0,
76  IterativeSolverResult *res = 0);
77 // PCG solver for multiple right-hand sides
78 
79 template<class MT>
80 int PCG (TVector<MT> (*Mv_clbk)(const TVector<MT> &v,
81  void *context), void *context, const TVector<MT> &b, TVector<MT> &x,
82  double &tol, TPreconditioner<MT> *precon = 0,
83  int maxit = 0);
84 
85 template<class MT>
86 int BiPCG (const TMatrix<MT> &A, const TVector<MT> &b, TVector<MT> &x,
87  double &tol, TPreconditioner<MT> *precon = 0, int maxit = 0);
88 // BiPCG (preconditioned bi-conjugate gradient) linear solver
89 // Solves Ax = b for given tolerance and preconditioner
90 
91 template<class MT>
92 int BiCGSTAB (const TMatrix<MT> &A, const TVector<MT> &b,
93  TVector<MT> &x, double &tol, TPreconditioner<MT> *precon = 0,
94  int maxit = 0);
95 
96 template<class MT>
97 void BiCGSTAB (const TMatrix<MT> &A, const TVector<MT> *b,
98  TVector<MT> *x, int nrhs, double tol, int maxit = 0,
99  TPreconditioner<MT> *precon = 0, IterativeSolverResult *res = 0);
100 // BiCGSTAB solver for multiple right-hand sides
101 
102 template<class MT>
103 int BiCGSTAB (TVector<MT> (*Mv_clbk)(const TVector<MT> &v,
104  void *context), void *context, const TVector<MT> &b, TVector<MT> &x,
105  double &tol, TPreconditioner<MT> *precon = 0,
106  int maxit = 0);
107 
108 template<class MT>
109 int BiCGSTAB (void (*MVM)(TVector<MT> &), const TVector<MT> &b,
110  TVector<MT> &x, double &tol, TPreconditioner<MT> *precon = 0,
111  int maxit = 0);
112 
113 template<class MT>
114 int GMRES (const TMatrix<MT> &A, const TVector<MT> &b, TVector<MT> &x,
115  double &tol, TPreconditioner<MT> *precon = 0, int restart = 10,
116  int maxit = 0, void (*clbk)(void*) = 0);
117 // GMRES (generalised minimum residual) linear solver
118 // Solves Ax = b for given tolerance and preconditioner
119 
139 template<class MT>
140 int GMRES (TVector<MT> (*Mv_clbk)(const TVector<MT> &v, void *context),
141  void *context, const TVector<MT> &b, TVector<MT> &x, double tol,
142  TPreconditioner<MT> *precon = 0, int restart = 10, int maxit = 0,
143  int *iter = 0, double *res = 0);
144 
145 template<class MT>
146 std::ostream &operator<< (std::ostream &os, const TMatrix<MT> &mat);
147 
148 // ==========================================================================
149 // class TMatrix
150 // ==========================================================================
161 template<class MT> class TMatrix {
162 public:
163  enum RC { ROW, COL }; // row/column identifier
164 
168  TMatrix ()
169  { rows = cols = 0; }
170 
178  TMatrix (int nrows, int ncols)
179  { rows = nrows, cols = ncols; }
180 
185  TMatrix (const TMatrix<MT> &m);
186 
191  virtual ~TMatrix() {}
192 
199  int Dim (RC rc) const { return (rc==ROW ? rows : cols); }
200 
206  int nRows () const { return rows; }
207 
213  int nCols () const { return cols; }
214 
219  bool isSparse() const { return StorageType() != MATRIX_DENSE; }
220 
223  bool isFull() const { return StorageType() == MATRIX_DENSE; }
224 
229  virtual MatrixStorage StorageType() const = 0;
230 
240  virtual void New (int nrows, int ncols);
241 
251  virtual MT Get (int r, int c) const = 0;
252 
260  inline MT operator () (int r, int c) const { return Get (r,c); }
261 
270  virtual TVector<MT> Row (int r) const = 0;
271 
283  virtual void SetRow (int r, const TVector<MT> &row)
284  { ERROR_UNDEF; }
285 
301  virtual int SparseRow (int r, idxtype *colidx, MT *val) const = 0;
302 
311  virtual TVector<MT> Col (int c) const = 0;
312 
321  virtual TVector<MT> Diag () const;
322 
334  virtual TVector<MT> ColNorm () const;
335 
336  virtual void ColScale (const TVector<MT> &scale) = 0;
337  // scales the columns with 'scale'
338 
339  virtual void RowScale (const TVector<MT> &scale) = 0;
340  // scales the rows with 'scale'
341 
342  virtual void Unlink () = 0;
343  // removes the matrix' link to its data block and deletes the data
344  // block, if necessary
345 
346  virtual void Ax (const TVector<MT> &x, TVector<MT> &b) const;
347  // returns the result of Ax in b; this is an alternative to
348  // operator* which avoids local vector creation and copy
349 
350  inline TVector<MT> operator* (const TVector<MT> &x) const
351  { TVector<MT> b; Ax (x, b); return b; }
352  // matrix x vector multiplication
353 
354  virtual void ATx (const TVector<MT> &x, TVector<MT> &b) const;
355  // returns the result of A^T x in b
356 
357  inline TVector<MT> ATx (const TVector<MT> &x) const
358  { TVector<MT> b; ATx (x, b); return b; }
359 
360  virtual void Transpone ()
361  { ERROR_UNDEF; }
362  // Replace matrix with its transpose
363 
364  virtual MT RowMult (int r, MT *x) const
365  { ERROR_UNDEF; return 0; }
366  // inner product of row r with full vector given by x
367  // where dimension(x) >= ncols
368 
376  void Export (std::ostream &os) const;
377 
378  void Print (std::ostream &os = std::cout, int n = 80) const;
379  // Pretty-print matrix in dense format to os, using a maximum of
380  // n characters per line
381 
382  void PrintNzeroGraph (char *fname);
383  // Generates a PPM bitmap file 'fname' with the image of the nonzero
384  // elements of the matrix.
385 
386  // **** friends ****
387 
389  friend TSymMatrix<MT> ATA<> (const TMatrix<MT> &A);
390  friend TSymMatrix<MT> AAT<> (const TMatrix<MT> &A);
391  friend TVector<MT>ATA_diag<> (const TMatrix<MT> &A);
392 
393  friend int PCG<> (const TMatrix<MT> &A, const TVector<MT> &b,
394  TVector<MT> &x, double &tol, TPreconditioner<MT> *precon,
395  int maxit);
396 
397  friend void PCG<> (const TMatrix<MT> &A, const TVector<MT> *b,
398  TVector<MT> *x, int nrhs, double tol, int maxit,
399  TPreconditioner<MT> *precon, IterativeSolverResult *res);
400 
401  friend int BiCGSTAB<> (const TMatrix<MT> &A,
402  const TVector<MT> &b, TVector<MT> &x, double &tol,
403  TPreconditioner<MT> *precon, int maxit);
404  // biconjugate gradient stabilised method.
405 
406  friend void BiCGSTAB<> (const TMatrix<MT> &A,
407  const TVector<MT> *b, TVector<MT> *x, int nrhs, double tol, int maxit,
408  TPreconditioner<MT> *precon, IterativeSolverResult *res);
409 
410  // Explicit linear solvers
411 
412  virtual int pcg (const TVector<MT> &b, TVector<MT> &x,
413  double &tol, TPreconditioner<MT> *precon = 0, int maxit = 0)
414  const;
415 
416  virtual void pcg (const TVector<MT> *b, TVector<MT> *x,
417  int nrhs, double tol, int maxit = 0,
418  TPreconditioner<MT> *precon = 0, IterativeSolverResult *res = 0)
419  const;
420 
421  virtual int bicgstab (const TVector<MT> &b, TVector<MT> &x,
422  double &tol, TPreconditioner<MT> *precon = 0, int maxit = 0)
423  const;
424 
425  virtual void bicgstab (const TVector<MT> *b, TVector<MT> *x,
426  int nrhs, double tol, int maxit = 0,
427  TPreconditioner<MT> *precon = 0, IterativeSolverResult *res = 0)
428  const;
429 
430  friend std::ostream &operator<< <> (std::ostream &os,
431  const TMatrix<MT> &mat);
432  // stream output
433 
434 protected:
435  int rows, cols; // number of rows and columns of the matrix
436 };
437 
438 // ==========================================================================
439 // typedefs for specific instances of `TMatrix'
440 
441 typedef TMatrix<double> RMatrix; // 'double real'
442 typedef TMatrix<float> FMatrix; // 'single real'
443 typedef TMatrix<std::complex<double> > CMatrix; // 'complex'
444 typedef TMatrix<std::complex<float> > SCMatrix; // 'single complex'
445 typedef TMatrix<int> IMatrix; // 'integer'
446 
447 
448 // ==========================================================================
449 // ==========================================================================
450 // Member definitions
451 
452 template<class MT>
453 TMatrix<MT>::TMatrix (const TMatrix<MT> &m)
454 {
455  rows = m.rows;
456  cols = m.cols;
457 }
458 
459 // --------------------------------------------------------------------------
460 
461 template<class MT>
462 void TMatrix<MT>::New (int nrows, int ncols)
463 {
464  rows = nrows;
465  cols = ncols;
466 }
467 
468 // --------------------------------------------------------------------------
469 
470 template<class MT>
471 TVector<MT> TMatrix<MT>::Diag () const
472 {
473  int i, n = (rows < cols ? rows : cols);
474  TVector<MT> diag(n);
475  for (i = 0; i < n; i++)
476  diag[i] = Get(i,i);
477  return diag;
478 }
479 
480 // --------------------------------------------------------------------------
481 
482 template<class MT>
483 TVector<MT> TMatrix<MT>::ColNorm () const
484 {
485  int j;
486  TVector<MT> cn(cols);
487  for (j = 0; j < cols; j++) {
488  cn[j] = (MT)l2norm (Col(j));
489  }
490  return cn;
491 }
492 
493 // --------------------------------------------------------------------------
494 
495 template<class MT>
496 void TMatrix<MT>::Ax (const TVector<MT> &x, TVector<MT> &b) const
497 {
498  dASSERT (cols == x.Dim(), "Argument 1: vector has wrong size");
499  if (b.Dim() != rows) b.New(rows);
500  for (int i = 0; i < rows; i++) b[i] = Row(i) & x;
501 }
502 
503 // --------------------------------------------------------------------------
504 
505 template<class MT>
506 void TMatrix<MT>::ATx (const TVector<MT> &x, TVector<MT> &b) const
507 {
508  dASSERT (rows == x.Dim(), "Argument 1: vector has wrong size");
509  if (b.Dim() != cols) b.New(cols);
510  for (int i = 0; i < cols; i++) b[i] = Col(i) & x;
511 }
512 
513 // --------------------------------------------------------------------------
514 
515 template<class MT>
516 TSymMatrix<MT> ATA (const TMatrix<MT> &A)
517 {
518  int i, j, nr = A.nRows(), nc = A.nCols();
519  TSymMatrix<MT> ata(nc);
520  for (i = 0; i < nc; i++) {
521  TVector<MT> col = A.Col(i);
522  for (j = 0; j <= i; j++)
523  ata(i,j) = col & A.Col(j);
524  }
525  return ata;
526 }
527 
528 // --------------------------------------------------------------------------
529 
530 template<>
531 inline int TMatrix<float>::pcg (const FVector &b, FVector &x,
532  double &tol, TPreconditioner<float> *precon, int maxit) const
533 {
534  return PCG (*this, b, x, tol, precon, maxit);
535 }
536 
537 template<>
538 inline int TMatrix<double>::pcg (const RVector &b, RVector &x,
539  double &tol, TPreconditioner<double> *precon, int maxit) const
540 {
541  return PCG (*this, b, x, tol, precon, maxit);
542 }
543 
544 template<class MT>
545 int TMatrix<MT>::pcg (const TVector<MT> &b, TVector<MT> &x,
546  double &tol, TPreconditioner<MT> *precon, int maxit) const
547 {
548  ERROR_UNDEF;
549  return 0;
550 }
551 
552 // --------------------------------------------------------------------------
553 
554 template<>
555 inline void TMatrix<float>::pcg (const FVector *b, FVector *x, int nrhs,
556  double tol, int maxit, TPreconditioner<float> *precon,
557  IterativeSolverResult *res) const
558 {
559  PCG (*this, b, x, nrhs, tol, maxit, precon, res);
560 }
561 
562 template<>
563 inline void TMatrix<double>::pcg (const RVector *b, RVector *x, int nrhs,
564  double tol, int maxit, TPreconditioner<double> *precon,
565  IterativeSolverResult *res) const
566 {
567  PCG (*this, b, x, nrhs, tol, maxit, precon, res);
568 }
569 
570 template<class MT>
571 void TMatrix<MT>::pcg (const TVector<MT> *b, TVector<MT> *x, int nrhs,
572  double tol, int maxit, TPreconditioner<MT> *precon,
573  IterativeSolverResult *res) const
574 {
575  ERROR_UNDEF;
576 }
577 
578 // --------------------------------------------------------------------------
579 
580 template<>
581 inline int TMatrix<float>::bicgstab (const FVector &b, FVector &x,
582  double &tol, TPreconditioner<float> *precon, int maxit) const
583 {
584  return BiCGSTAB (*this, b, x, tol, precon, maxit);
585 }
586 
587 template<>
588 inline int TMatrix<double>::bicgstab (const RVector &b, RVector &x,
589  double &tol, TPreconditioner<double> *precon, int maxit) const
590 {
591  return BiCGSTAB (*this, b, x, tol, precon, maxit);
592 }
593 
594 template<class MT>
595 int TMatrix<MT>::bicgstab (const TVector<MT> &b, TVector<MT> &x,
596  double &tol, TPreconditioner<MT> *precon, int maxit) const
597 {
598  ERROR_UNDEF;
599  return 0;
600 }
601 
602 // --------------------------------------------------------------------------
603 
604 template<>
605 inline void TMatrix<float>::bicgstab (const FVector *b, FVector *x, int nrhs,
606  double tol, int maxit,TPreconditioner<float> *precon,
607  IterativeSolverResult *res) const
608 {
609  BiCGSTAB (*this, b, x, nrhs, tol, maxit, precon, res);
610 }
611 
612 template<>
613 inline void TMatrix<double>::bicgstab (const RVector *b, RVector *x, int nrhs,
614  double tol, int maxit, TPreconditioner<double> *precon,
615  IterativeSolverResult *res) const
616 {
617  BiCGSTAB (*this, b, x, nrhs, tol, maxit, precon, res);
618 }
619 
620 template<class MT>
621 void TMatrix<MT>::bicgstab (const TVector<MT> *b, TVector<MT> *x, int nrhs,
622  double tol, int maxit, TPreconditioner<MT> *precon,
623  IterativeSolverResult *res) const
624 {
625  ERROR_UNDEF;
626 }
627 
628 // --------------------------------------------------------------------------
629 
630 template<class MT>
631 void TMatrix<MT>::Export (std::ostream &os) const
632 {
633  int r, c;
634  for (r = 0; r < rows; r++) {
635  for (c = 0; c < cols; c++) {
636  os << Get(r,c);
637  os << (c == cols-1 ? '\n' : ' ');
638  }
639  }
640 }
641 
642 // --------------------------------------------------------------------------
643 
644 template<class MT>
645 void TMatrix<MT>::Print (std::ostream &os, int n) const
646 {
647  int r, c, nc, minc, maxc, maxlen = 0;
648 
649  maxlen = 11;
650  nc = n / (maxlen+1);
651  if (cols <= nc) {
652  for (r = 0; r < rows; r++) {
653  for (c = 0; c < cols; c++) {
654  os.width(maxlen);
655  os.precision(5);
656  os << Get(r,c) << ' ';
657  }
658  os << std::endl;
659  }
660  } else {
661  minc = 0;
662  while (minc < cols) {
663  maxc = minc+nc;
664  if (maxc > cols) maxc = cols;
665  os << "*** Columns " << (minc+1) << " through " << maxc
666  << " ***\n";
667  for (r = 0; r < rows; r++) {
668  for (c = minc; c < maxc; c++) {
669  os.width(maxlen);
670  os.precision(5);
671  os << Get(r,c) << ' ';
672  }
673  os << std::endl;
674  }
675  minc = maxc;
676  }
677  }
678 }
679 
680 // --------------------------------------------------------------------------
681 
682 template<class MT>
683 void TMatrix<MT>::PrintNzeroGraph (char *fname)
684 {
685  int i, r, c;
686  std::ofstream ofs (fname);
687  ofs << "P1" << std::endl;
688  ofs << "# CREATOR: TOAST TMatrix::PrintNzeroGraph" << std::endl;
689  ofs << cols << ' ' << rows << std::endl;
690  for (i = r = 0; r < rows; r++) {
691  for (c = 0; c < cols; c++) {
692  ofs << (Get(r,c) != (MT)0 ? '1' : '0');
693  ofs << (++i % 35 ? ' ' : '\n');
694  }
695  }
696 }
697 
698 // --------------------------------------------------------------------------
699 
700 #endif // !__MATRIX_H
TMatrix::SetRow
virtual void SetRow(int r, const TVector< MT > &row)
Substitute a row of the matrix.
Definition: matrix.h:283
TVector::New
void New(int dim)
Resize the vector.
Definition: vector.h:338
MATRIX_SYMCOMPROW
symmetric compressed row storage (sparse)
Definition: matrix.h:27
TMatrix::isFull
bool isFull() const
Return dense storage flag.
Definition: matrix.h:223
MATRIX_COMPCOL
compressed column storage (sparse)
Definition: matrix.h:26
TMatrix::ATA
friend TSymMatrix< MT > ATA(const TMatrix< MT > &A)
Return transp(*this) * *this as a symmetric matrix.
Definition: matrix.h:516
TSymMatrix
Definition: dnsmatrix.h:17
TMatrix::TMatrix
TMatrix()
Create a matrix of size 0 x 0.
Definition: matrix.h:168
MATRIX_COMPROW
compressed row storage (sparse)
Definition: matrix.h:25
TPreconditioner
Definition: gsmatrix.h:48
TMatrix::SparseRow
virtual int SparseRow(int r, idxtype *colidx, MT *val) const =0
Returns a row of the matrix in sparse format.
TMatrix::New
virtual void New(int nrows, int ncols)
Resize and reset the matrix.
Definition: matrix.h:462
TMatrix::Col
virtual TVector< MT > Col(int c) const =0
Returns a vector containing a copy of column 'c'.
MatrixStorage
MatrixStorage
Definition: matrix.h:20
TMatrix::TMatrix
TMatrix(int nrows, int ncols)
Create a matrix of logical size nrows x ncols.
Definition: matrix.h:178
TMatrix::operator()
MT operator()(int r, int c) const
Matrix element access (read only)
Definition: matrix.h:260
MATRIX_DENSE
dense matrix storage
Definition: matrix.h:21
TVector::Dim
int Dim() const
Returns the size of the vector.
Definition: vector.h:295
TMatrix::ColNorm
virtual TVector< MT > ColNorm() const
Returns vector of column norms.
Definition: matrix.h:483
TMatrix
Virtual base class for all matrix types (dense and sparse)
Definition: matrix.h:43
TMatrix::nCols
int nCols() const
Return number of columns of the matrix.
Definition: matrix.h:213
TMatrix::Row
virtual TVector< MT > Row(int r) const =0
Returns a vector containing a copy of row `r'.
TMatrix::Get
virtual MT Get(int r, int c) const =0
Retrieve a matrix element.
TMatrix::~TMatrix
virtual ~TMatrix()
Destroy the matrix.
Definition: matrix.h:191
TMatrix::Export
void Export(std::ostream &os) const
Write matrix to ASCII stream.
Definition: matrix.h:631
TMatrix::nRows
int nRows() const
Return number of rows of the matrix.
Definition: matrix.h:206
TMatrix::isSparse
bool isSparse() const
Return sparse storage flag.
Definition: matrix.h:219
IterativeSolverResult
Definition: matrix.h:34
TMatrix::Dim
int Dim(RC rc) const
Return a matrix dimension.
Definition: matrix.h:199
MATRIX_DIAG
store diagonal matrix as vector
Definition: matrix.h:23
TMatrix::StorageType
virtual MatrixStorage StorageType() const =0
Matrix storage class.
MATRIX_COORD
coordinate storage (sparse)
Definition: matrix.h:24
MATRIX_SYM
store lower triangle + diagonal of symmetric matrix
Definition: matrix.h:22
TMatrix::Diag
virtual TVector< MT > Diag() const
Returns the matrix diagonal as a vector.
Definition: matrix.h:471
Generated on Tue Sep 2 2014 17:14:26 for Toast++ by   doxygen 1.8.6