Matrix.h 21.2 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17 18 19 20 21 22 23
// Matrix.h
//
////////////////////////////////////////////////////////////////////////////////

#ifndef MATRIX_H
#define MATRIX_H

#define MATRIX_DEF_BLOCK_SIZE 64

24
#include "Context.h"
Francois Gygi committed
25 26 27

#include <valarray>
#include <complex>
28
#include <cstring> // memcpy
Francois Gygi committed
29 30 31 32 33 34

class ComplexMatrix;

class DoubleMatrix
{
  private:
35

Francois Gygi committed
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
    Context ctxt_;
    int ictxt_;
    int lld_;  // leading dimension of local matrix
    int m_, n_;   // size of global matrix
    int mb_, nb_; // size of blocks
    int mloc_, nloc_;   // size of local array
    int size_; // size of local array;
    int nprow_, npcol_; // number of process rows and cols in context
    int myrow_, mycol_; // position of my process in the process grid
    int mblocks_, nblocks_; // number of local blocks
    int desc_[9];
    bool active_;
    bool m_incomplete_, n_incomplete_; // this process has an incomplete block
    bool reference_; // object was created using the copy constructor
    double* val;

  public:
53

Francois Gygi committed
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
    double* valptr(int i=0) { return &val[i]; }
    const double* cvalptr(int i=0) const { return &val[i]; }
    double& operator[] (int i) { return val[i]; }
    const double& operator[] (int i) const { return val[i]; }
    const Context& context(void) const { return ctxt_; }
    int ictxt(void) const { return ictxt_; }
    int lld(void) const { return lld_; }
    int m(void) const { return m_; } // size of global matrix
    int n(void) const { return n_; } // size of global matrix
    int mb(void) const { return mb_; } // size of blocks
    int nb(void) const { return nb_; } // size of blocks
    int mloc(void) const { return mloc_; } // size of local array
    int nloc(void) const { return nloc_; } // size of local array
    int size(void) const { return size_; } // size of local array
    int localsize(void) const { return mloc_*nloc_; } // local size of val
    double memsize(void) const { return (double)m_ * (double)n_ * sizeof(double); }
    double localmemsize(void) const
    { return (double) mloc_ * (double) nloc_ * sizeof(double); }
    const int* desc(void) const { return &desc_[0]; }
73

Francois Gygi committed
74 75 76 77 78 79 80 81 82 83 84 85 86
    // local block size of block (l,m)
    int mbs(int l) const
    {
      // if l is the last block and this process holds an incomplete
      // block, then the size is m_%mb_. Otherwise, it is mb_.
      return ( (l==mblocks_-1) && ( m_incomplete_ ) ) ? m_%mb_ : mb_;
    }
    int nbs(int m) const
    {
      // if m is the last block and this process holds an incomplete
      // block, then the size is n_%nb_. Otherwise, it is nb_.
      return ( (m==nblocks_-1) && ( n_incomplete_ ) ) ? n_%nb_ : nb_;
    }
87

Francois Gygi committed
88
    // number of local blocks
89
    int mblocks(void) const { return mblocks_; }
Francois Gygi committed
90
    int nblocks(void) const { return nblocks_; }
91

Francois Gygi committed
92 93 94
    // functions to compute local indices from global indices

    // index of blocks: element (i,j) is at position (x,y)
95
    // in local block (l,m) of process (pr,pc)
Francois Gygi committed
96 97 98
    int l(int i) const { return i/(nprow_*mb_); }
    int x(int i) const { return i % mb_; }
    int pr(int i) const { return (i/mb_) % nprow_; }
99

Francois Gygi committed
100 101 102
    int m(int j) const { return j/(npcol_*nb_); }
    int y(int j) const { return j % nb_; }
    int pc(int j) const { return (j/nb_) % npcol_; }
103

Francois Gygi committed
104 105
    // global indices:
    // (i,j) is the global index of element (x,y) of block (l,m)
106
    int i(int l, int x) const { return (l * nprow_ + myrow_) * mb_ + x; }
Francois Gygi committed
107
    int j(int m, int y) const { return (m * npcol_ + mycol_) * nb_ + y; }
108 109 110 111 112

    int iglobal(int ilocal) const
    { return mb_*(nprow_*(ilocal/mb_)+myrow_)+ilocal%mb_; }
    int jglobal(int jlocal) const
    { return nb_*(npcol_*(jlocal/nb_)+mycol_)+jlocal%nb_; }
113

114 115 116 117 118 119
    int iglobal(int iprow, int ilocal) const
    { return mb_*(nprow_*(ilocal/mb_)+iprow)+ilocal%mb_; }
    int jglobal(int jpcol, int jlocal) const
    { return nb_*(npcol_*(jlocal/nb_)+jpcol)+jlocal%nb_; }


Francois Gygi committed
120 121 122 123 124
    // store element a(ii,jj) (where ii,jj are global indices)
    // in array val:
    //        int iii = l(ii) * mb_ + x(ii);
    //        int jjj = m(jj) * nb_ + y(jj);
    //        val[iii+mloc_*jjj] = a_ij;
125

Francois Gygi committed
126 127
    // active flag: the matrix has elements on this process
    bool active(void) const { return active_; }
128 129

    void init_size(int m, int n, int mb = MATRIX_DEF_BLOCK_SIZE,
Francois Gygi committed
130
      int nb = MATRIX_DEF_BLOCK_SIZE);
131 132

    void resize(int m, int n, int mb = MATRIX_DEF_BLOCK_SIZE,
Francois Gygi committed
133 134
      int nb = MATRIX_DEF_BLOCK_SIZE)
    {
135 136 137
      const int old_size = size_;
      init_size(m,n,mb,nb);
      if ( size_ == old_size ) return;
Francois Gygi committed
138 139
      delete[] val;
      val = new double[size_];
140
      clear();
141 142
    }

143
    void print(std::ostream& os) const;
144

Francois Gygi committed
145
    explicit DoubleMatrix(const Context& ctxt) : ctxt_(ctxt),
146 147
        m_(0), n_(0), mb_(0), nb_(0), mloc_(0), nloc_(0),
        size_(0), reference_(false), val(0) {}
148

Francois Gygi committed
149 150
    // Construct a DoubleMatrix of dimensions m,n
    explicit DoubleMatrix(const Context& ctxt, int m, int n,
151
        int mb=MATRIX_DEF_BLOCK_SIZE,
152
        int nb=MATRIX_DEF_BLOCK_SIZE) : ctxt_(ctxt),
153
        m_(0), n_(0), mb_(0), nb_(0), size_(0), reference_(false), val(0)
Francois Gygi committed
154 155 156
    {
      resize(m,n,mb,nb);
    }
157

Francois Gygi committed
158 159 160 161 162 163 164 165
    // copy constructor: create a separate copy of rhs
    explicit DoubleMatrix(const DoubleMatrix& rhs) : ctxt_(rhs.context()),
      reference_(false)
    {
      init_size(rhs.m(),rhs.n(),rhs.mb(),rhs.nb());
      val = new double[size_];
      memcpy(val, rhs.val, size_*sizeof(double));
    }
166

Francois Gygi committed
167 168 169
    // reference constructor create a proxy for a ComplexMatrix rhs
    explicit DoubleMatrix(ComplexMatrix& rhs);
    explicit DoubleMatrix(const ComplexMatrix& rhs);
170 171 172

    ~DoubleMatrix(void)
    {
Francois Gygi committed
173 174
      if ( !reference_ ) delete[] val;
    }
175

Francois Gygi committed
176
    DoubleMatrix& operator=(const DoubleMatrix& a);
177

Francois Gygi committed
178 179
    DoubleMatrix& operator+=(const DoubleMatrix& a);
    DoubleMatrix& operator-=(const DoubleMatrix& a);
180

Francois Gygi committed
181
    DoubleMatrix& operator*=(double a);
182

Francois Gygi committed
183 184 185 186 187 188 189
    void matgather(double *a, int lda) const;

    void initdiag(const double* const a); // initialize diagonal from a[]
    void init(const double* const a, int lda);
    void set(char uplo, double x);
    void identity(void);
    void clear(void);
190

Francois Gygi committed
191 192 193
    double dot(const DoubleMatrix &a) const;
    void axpy(double alpha, const DoubleMatrix &a);
    void scal(double alpha);
194

Francois Gygi committed
195 196 197 198 199
    double nrm2(void) const;
    double asum(void) const;
    double amax(void) const;
    double trace(void) const;
    void symmetrize(char uplo);
200

Francois Gygi committed
201 202 203 204 205 206 207 208
    // rank-1 update: *this += alpha * x(kx) * (y(ky))^T
    // where x(kx) is row kx of x, and y(ky) is row ky of y
    void ger(double alpha, const DoubleMatrix& x, int kx,
                           const DoubleMatrix& y, int ky);

    // symmetric rank-1 update
    void syr(char uplo, double alpha, const DoubleMatrix& x,
             int ix, char rowcol);
209

Francois Gygi committed
210 211
    // get submatrix A(ia:ia+m,ja:ja+n) of A
    void getsub(const DoubleMatrix& a,int m,int n,int ia,int ja);
212

213 214 215 216 217
    // get submatrix A(ia:ia+m,ja:ja+n) of A and store in
    // this(idest:idest+m,jdest:jdest+n)
    void getsub(const DoubleMatrix& a,int m,int n,int isrc,int jsrc,
                int idest, int jdest);

Francois Gygi committed
218 219 220 221 222
    // matrix * matrix
    // this = alpha*op(A)*op(B)+beta*this
    void gemm(char transa, char transb,
         double alpha, const DoubleMatrix& a,
         const DoubleMatrix& b, double beta);
223

Francois Gygi committed
224 225 226 227 228
    // symmetric_matrix * matrix
    // *this = alpha * A * B + beta * this
    void symm(char side, char uplo,
         double alpha, const DoubleMatrix& a,
         const DoubleMatrix& b, double beta);
229

Francois Gygi committed
230 231 232 233 234 235 236 237 238 239 240
    // symmetric rank k update
    // this = beta * this + alpha * A * A^T  (trans=='n')
    // this = beta * this + alpha * A^T * A  (trans=='t')
    void syrk(char uplo, char trans,
         double alpha, const DoubleMatrix& a, double beta);

    // matrix transpose
    // this = alpha * transpose(A) + beta * this
    void transpose(double alpha, const DoubleMatrix& a, double beta);
    void transpose(const DoubleMatrix& a);

241
    void trmm(char side, char uplo, char trans, char diag,
Francois Gygi committed
242
              double alpha, const DoubleMatrix& a);
243

244
    // solve triangular system
245
    void trsm(char side, char uplo, char trans, char diag,
Francois Gygi committed
246 247
              double alpha, const DoubleMatrix& a);
    void trtrs(char uplo, char trans, char diag, DoubleMatrix& b) const;
248

Francois Gygi committed
249 250 251 252
    // Cholesky decomposition of a symmetric matrix
    void potrf(char uplo);
    // Inverse of a symmetric matrix from Cholesky factor
    void potri(char uplo);
253

254 255
    // Polar decomposition, tolerance ||I-X^T*X||<tol or iter<maxiter
    void polar(double tol, int maxiter);
Francois Gygi committed
256

257
    // LU decomposition
258
    void lu(std::valarray<int>& ipiv);
259

260 261
    // compute inverse of a square matrix
    void inverse(void);
262 263
    // compute inverse and determinant of a square matrix
    double inverse_det(void);
264 265 266
    // compute determinant of a square matrix in LU form
    double det_from_lu(std::valarray<int> ipiv);
    // compute inverse of a square matrix in LU form
267
    void inverse_from_lu(std::valarray<int>& ipiv);
268

Francois Gygi committed
269 270
    // Inverse of triangular matrix
    void trtri(char uplo,char diag);
271

Francois Gygi committed
272 273
    double pocon(char) const;
    void sygst(int, char, const DoubleMatrix&);
274

Francois Gygi committed
275
    // compute eigenvalues and eigenvectors of symmetric matrix *this
276
    void syev(char uplo, std::valarray<double>& w, DoubleMatrix& z);
277

Francois Gygi committed
278
    // compute eigenvalues (only) of symmetric matrix *this
279
    void syev(char uplo, std::valarray<double>& w);
280 281 282

    // compute eigenvalues and eigenvectors of symmetric matrix *this
    // using the divide and conquer method of Tisseur and Dongarra
283
    void syevd(char uplo, std::valarray<double>& w, DoubleMatrix& z);
284

285 286
    // compute eigenvalues (only) of symmetric matrix *this
    // using the divide and conquer method of Tisseur and Dongarra
287
    void syevd(char uplo, std::valarray<double>& w);
288 289 290

    // compute eigenvalues and eigenvectors of symmetric matrix *this
    // using the expert driver
291
    void syevx(char uplo, std::valarray<double>& w, DoubleMatrix& z,
292
       double abstol);
293

294
    // permute the coeff of the matrix *this
295
    void lapiv(char direc, char rowcol, const int *ipiv);
296 297
    // signature of a permutation returned by lu
    int signature(std::valarray<int> ipiv);
298

299 300
    // compute eigenvalues (only) of symmetric matrix *this
    // using the divide and conquer method of Tisseur and Dongarra
301
    //void syevx(char uplo, std::valarray<double>& w);
Francois Gygi committed
302
};
303
std::ostream& operator << ( std::ostream& os, const DoubleMatrix& a );
Francois Gygi committed
304 305 306 307

class ComplexMatrix
{
  private:
308

Francois Gygi committed
309 310 311 312 313 314 315 316 317 318 319 320 321 322
    Context ctxt_;
    int ictxt_;
    int lld_;  // leading dimension of local matrix
    int m_, n_;   // size of global matrix
    int mb_, nb_; // size of blocks
    int mloc_, nloc_;   // size of local array
    int size_; // size of local array;
    int nprow_, npcol_; // number of process rows and cols in  context
    int myrow_, mycol_; // position of my process in the process grid
    int mblocks_, nblocks_; // number of local blocks
    int desc_[9];
    bool active_;
    bool m_incomplete_, n_incomplete_; // this process has an incomplete block
    bool reference_; // object was created using the copy constructor
323
    std::complex<double>* val;
Francois Gygi committed
324 325

  public:
326

327 328 329 330
    std::complex<double>* valptr(int i=0) { return &val[i]; }
    const std::complex<double>* cvalptr(int i=0) const { return &val[i]; }
    std::complex<double>& operator[] (int i) { return val[i]; }
    const std::complex<double>& operator[] (int i) const { return val[i]; }
Francois Gygi committed
331 332 333 334 335 336 337 338 339 340 341
    const Context& context(void) const { return ctxt_; }
    int ictxt(void) const { return ictxt_; }
    int lld(void) const { return lld_; }
    int m(void) const { return m_; } // size of global matrix
    int n(void) const { return n_; } // size of global matrix
    int mb(void) const { return mb_; } // size of blocks
    int nb(void) const { return nb_; } // size of blocks
    int mloc(void) const { return mloc_; } // size of local array
    int nloc(void) const { return nloc_; } // size of local array
    int size(void) const { return size_; } // size of local array
    int localsize(void) const { return mloc_*nloc_; } // local size of val
342
    double memsize(void) const
343
    { return (double) m_ * (double) n_ * sizeof(std::complex<double>); }
Francois Gygi committed
344
    double localmemsize(void) const
345
    { return (double) mloc_ * (double) nloc_ * sizeof(std::complex<double>); }
Francois Gygi committed
346
    const int* desc(void) const { return &desc_[0]; }
347

Francois Gygi committed
348 349 350 351 352 353 354 355 356 357 358 359 360
    // local block size of block (l,m)
    int mbs(int l) const
    {
      // if l is the last block and this process holds an incomplete
      // block, then the size is m_%mb_. Otherwise, it is mb_.
      return ( (l==mblocks_-1) && ( m_incomplete_ ) ) ? m_%mb_ : mb_;
    }
    int nbs(int m) const
    {
      // if m is the last block and this process holds an incomplete
      // block, then the size is n_%nb_. Otherwise, it is nb_.
      return ( (m==nblocks_-1) && ( n_incomplete_ ) ) ? n_%nb_ : nb_;
    }
361

Francois Gygi committed
362
    // number of local blocks
363
    int mblocks(void) const { return mblocks_; }
Francois Gygi committed
364
    int nblocks(void) const { return nblocks_; }
365

Francois Gygi committed
366 367 368
    // functions to compute local indices from global indices

    // index of blocks: element (i,j) is at position (x,y)
369
    // in local block (l,m) of process (pr,pc)
Francois Gygi committed
370 371 372
    int l(int i) const { return i/(nprow_*mb_); }
    int x(int i) const { return i % mb_; }
    int pr(int i) const { return (i/mb_) % nprow_; }
373

Francois Gygi committed
374 375 376
    int m(int j) const { return j/(npcol_*nb_); }
    int y(int j) const { return j % nb_; }
    int pc(int j) const { return (j/nb_) % npcol_; }
377

Francois Gygi committed
378 379
    // global indices:
    // (i,j) is the global index of element (x,y) of block (l,m)
380
    int i(int l, int x) const { return (l * nprow_ + myrow_) * mb_ + x; }
Francois Gygi committed
381
    int j(int m, int y) const { return (m * npcol_ + mycol_) * nb_ + y; }
382

383 384 385 386 387
    int iglobal(int ilocal) const
    { return mb_*(nprow_*(ilocal/mb_)+myrow_)+ilocal%mb_; }
    int jglobal(int jlocal) const
    { return nb_*(npcol_*(jlocal/nb_)+mycol_)+jlocal%nb_; }

388 389 390 391 392
    int iglobal(int iprow, int ilocal) const
    { return mb_*(nprow_*(ilocal/mb_)+iprow)+ilocal%mb_; }
    int jglobal(int jpcol, int jlocal) const
    { return nb_*(npcol_*(jlocal/nb_)+jpcol)+jlocal%nb_; }

Francois Gygi committed
393 394 395 396 397
    // store element a(ii,jj) (where ii,jj are global indices)
    // in array val:
    //        int iii = l(ii) * mb_ + x(ii);
    //        int jjj = m(jj) * nb_ + y(jj);
    //        val[iii+mloc_*jjj] = a_ij;
398

Francois Gygi committed
399 400
    // active flag: the matrix has elements on this process
    bool active(void) const { return active_; }
401 402

    void init_size(int m, int n, int mb = MATRIX_DEF_BLOCK_SIZE,
Francois Gygi committed
403
      int nb = MATRIX_DEF_BLOCK_SIZE);
404 405

    void resize(int m, int n, int mb = MATRIX_DEF_BLOCK_SIZE,
Francois Gygi committed
406 407
      int nb = MATRIX_DEF_BLOCK_SIZE)
    {
408
      const int old_size = size_;
Francois Gygi committed
409
      init_size(m,n,mb,nb);
410
      if ( size_ == old_size ) return;
Francois Gygi committed
411
      delete[] val;
412
      val = new std::complex<double>[size_];
413
      clear();
414 415
    }

416
    void print(std::ostream& os) const;
417

Francois Gygi committed
418
    explicit ComplexMatrix(const Context& ctxt) : ctxt_(ctxt),
419 420
        m_(0), n_(0), mb_(0), nb_(0), mloc_(0), nloc_(0),
        size_(0), reference_(false), val(0) {}
421

Francois Gygi committed
422 423
    // Construct a ComplexMatrix of dimensions m,n
    explicit ComplexMatrix(const Context& ctxt, int m, int n,
424 425
        int mb=MATRIX_DEF_BLOCK_SIZE,
        int nb=MATRIX_DEF_BLOCK_SIZE) : ctxt_(ctxt),
426
        m_(0), n_(0), mb_(0), nb_(0), size_(0), reference_(false), val(0)
Francois Gygi committed
427 428 429
    {
      resize(m,n,mb,nb);
    }
430

Francois Gygi committed
431 432 433 434 435
    // copy constructor: create a separate copy of rhs
    explicit ComplexMatrix(const ComplexMatrix& rhs) : ctxt_(rhs.context()),
      reference_(false)
    {
      init_size(rhs.m(),rhs.n(),rhs.mb(),rhs.nb());
436 437
      val = new std::complex<double>[size_];
      memcpy(val, rhs.val, size_*sizeof(std::complex<double>));
Francois Gygi committed
438
    }
439

Francois Gygi committed
440 441 442
    // reference constructor: create a proxy for a DoubleMatrix rhs
    explicit ComplexMatrix(DoubleMatrix& rhs);
    explicit ComplexMatrix(const DoubleMatrix& rhs);
443

Francois Gygi committed
444
    ~ComplexMatrix(void)
445
    {
Francois Gygi committed
446 447
      if ( !reference_ ) delete[] val;
    }
448

Francois Gygi committed
449
    ComplexMatrix& operator=(const ComplexMatrix& a);
450

Francois Gygi committed
451 452
    ComplexMatrix& operator+=(const ComplexMatrix& a);
    ComplexMatrix& operator-=(const ComplexMatrix& a);
453

Francois Gygi committed
454
    ComplexMatrix& operator*=(double a);
455
    ComplexMatrix& operator*=(std::complex<double> a);
456

457
    void matgather(std::complex<double> *a, int lda) const;
Francois Gygi committed
458

459
    void initdiag(const std::complex<double>* const a); // initialize diagonal
460 461
    void init(const std::complex<double>* const a, int lda);
    void set(char uplo, std::complex<double> x);
Francois Gygi committed
462 463
    void identity(void);
    void clear(void);
464

465 466 467
    std::complex<double> dot(const ComplexMatrix &a) const;  // tr A^H * A
    std::complex<double> dotu(const ComplexMatrix &a) const; // tr A^T * A
    void axpy(std::complex<double> alpha, const ComplexMatrix &a);
Francois Gygi committed
468
    void axpy(double alpha, const ComplexMatrix &a);
469
    void scal(std::complex<double> alpha);
Francois Gygi committed
470
    void scal(double alpha);
471

Francois Gygi committed
472 473 474
    double nrm2(void) const;
    double asum(void) const;
    double amax(void) const;
475
    std::complex<double> trace(void) const;
Francois Gygi committed
476
    void symmetrize(char uplo);
477

Francois Gygi committed
478 479
    // rank-1 update: *this += alpha * x(kx) * (y(ky))^T
    // where x(kx) is row kx of x, and y(ky) is row ky of y
480
    void ger(std::complex<double> alpha, const ComplexMatrix& x,int kx,
Francois Gygi committed
481
                                    const ComplexMatrix& y,int ky);
482
    void geru(std::complex<double> alpha, const ComplexMatrix& x,int kx,
Francois Gygi committed
483
                                     const ComplexMatrix& y,int ky);
484
    void gerc(std::complex<double> alpha, const ComplexMatrix& x,int kx,
Francois Gygi committed
485 486 487
                                     const ComplexMatrix& y,int ky);

    // symmetric rank-1 update
488
    void her(char uplo, std::complex<double> alpha,
Francois Gygi committed
489
             const ComplexMatrix& x, int ix, char rowcol);
490

Francois Gygi committed
491 492
    // get submatrix A(ia:ia+m,ja:ja+n) of A
    void getsub(const ComplexMatrix& a, int m, int n, int ia, int ja);
493

494 495 496 497 498
    // get submatrix A(ia:ia+m,ja:ja+n) of A and store in
    // this(idest:idest+m,jdest:jdest+n)
    void getsub(const ComplexMatrix& a,int m,int n,int isrc,int jsrc,
                int idest, int jdest);

Francois Gygi committed
499 500 501
    // matrix * matrix
    // this = alpha*op(A)*op(B)+beta*this
    void gemm(char transa, char transb,
502 503
         std::complex<double> alpha, const ComplexMatrix& a,
         const ComplexMatrix& b, std::complex<double> beta);
504

Francois Gygi committed
505 506 507
    // hermitian_matrix * matrix
    // *this = alpha * A * B + beta * this
    void hemm(char side, char uplo,
508 509
         std::complex<double> alpha, const ComplexMatrix& a,
         const ComplexMatrix& b, std::complex<double> beta);
510

Francois Gygi committed
511 512 513
    // complex_symmetric_matrix * matrix
    // *this = alpha * A * B + beta * this
    void symm(char side, char uplo,
514 515
         std::complex<double> alpha, const ComplexMatrix& a,
         const ComplexMatrix& b, std::complex<double> beta);
516

Francois Gygi committed
517
    // hermitian rank k update
518 519
    // this = beta * this + alpha * A * A^H  (trans=='n')
    // this = beta * this + alpha * A^H * A  (trans=='c')
Francois Gygi committed
520
    void herk(char uplo, char trans,
521
         double alpha, const ComplexMatrix& a, double beta);
Francois Gygi committed
522 523 524

    // matrix transpose
    // this = alpha * transpose(A) + beta * this
525
    void transpose(std::complex<double> alpha, const ComplexMatrix& a,
526
                   std::complex<double> beta);
Francois Gygi committed
527 528
    void transpose(const ComplexMatrix& a);

529
    void trmm(char side, char uplo, char trans, char diag,
530
              std::complex<double> alpha, const ComplexMatrix& a);
531
    void trsm(char side, char uplo, char trans, char diag,
532
              std::complex<double> alpha, const ComplexMatrix& a);
Francois Gygi committed
533
    void trtrs(char uplo, char trans, char diag, ComplexMatrix& b) const;
534

Francois Gygi committed
535 536 537 538
    // Cholesky decomposition of a hermitian matrix
    void potrf(char uplo);
    // Inverse of a symmetric matrix from Cholesky factor
    void potri(char uplo);
539

540 541 542
    // Polar decomposition, tolerance ||I-X^H*X||<tol or iter<maxiter
    void polar(double tol, int maxiter);

543 544 545 546 547 548 549
    // LU decomposition
    void lu(std::valarray<int>& ipiv);

    // compute inverse of a square matrix
    void inverse(void);
    // compute inverse and determinant of a square matrix
    std::complex<double> inverse_det(void);
550 551 552
    // compute determinant of a square matrix in LU form
    std::complex<double> det_from_lu(std::valarray<int> ipiv);
    // compute inverse of a square matrix in LU form
553 554
    void inverse_from_lu(std::valarray<int>& ipiv);

Francois Gygi committed
555 556
    // Inverse of triangular matrix
    void trtri(char uplo,char diag);
557

Francois Gygi committed
558
    double pocon(char) const;
559

Francois Gygi committed
560
    // compute eigenvalues and eigenvectors of hermitian matrix *this
561
    void heev(char uplo, std::valarray<double>& w, ComplexMatrix& z);
Francois Gygi committed
562
    // compute eigenvalues (only) of hermitian matrix *this
563
    void heev(char uplo, std::valarray<double>& w);
Francois Gygi committed
564 565 566 567 568 569

    // eigenvalues/eigenvectors using the divide and conquer method
    // compute eigenvalues and eigenvectors of hermitian matrix *this
    void heevd(char uplo, std::valarray<double>& w, ComplexMatrix& z);
    // compute eigenvalues (only) of hermitian matrix *this
    void heevd(char uplo, std::valarray<double>& w);
570 571

    // permute the coeff of the matrix *this
572
    void lapiv(char direc, char rowcol, const int *ipiv);
573 574
    // signature of a permutation returned by lu
    int signature(std::valarray<int> ipiv);
Francois Gygi committed
575
};
576
std::ostream& operator << ( std::ostream& os, const ComplexMatrix& a );
Francois Gygi committed
577
#endif