Matrix.h

Go to the documentation of this file.
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure
00002  * calculations.
00003  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
00004  * 
00005  * This program is free software: you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation, either version 3 of the License, or
00008  * (at your option) any later version.
00009  * 
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  * 
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00017  * 
00018  * Primary academic reference:
00019  * Kohn−Sham Density Functional Theory Electronic Structure Calculations 
00020  * with Linearly Scaling Computational Time and Memory Usage,
00021  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
00022  * J. Chem. Theory Comput. 7, 340 (2011),
00023  * <http://dx.doi.org/10.1021/ct100611z>
00024  * 
00025  * For further information about Ergo, see <http://www.ergoscf.org>.
00026  */
00027 
00058 #ifndef MAT_MATRIX
00059 #define MAT_MATRIX
00060 
00061 #include <math.h>
00062 #include <cstdlib>
00063 #include <algorithm>
00064 
00065 #include "MatrixHierarchicBase.h"
00066 #include "matrix_proxy.h"
00067 #include "gblas.h"
00068 #include "sort.h"
00069 //#define MAT_NAIVE
00070 
00071 namespace mat{
00072   enum side {left, right};
00073   enum inchversion {unstable, stable};
00074   template<class Treal, class Telement>
00075     class Vector;
00076   template<typename Tperm>
00077     struct AccessMap;
00078 
00079   class SingletonForTimings {
00080   private:
00081     double accumulatedTime;
00082   public:
00083     void reset() { accumulatedTime = 0; }
00084     double getAccumulatedTime() { return accumulatedTime; }
00085     void addTime(double timeToAdd) { 
00086 #ifdef _OPENMP
00087       #pragma omp critical
00088 #endif
00089       {
00090         accumulatedTime += timeToAdd;
00091       }
00092     }
00093     static SingletonForTimings & instance() {
00094       static SingletonForTimings theInstance;
00095       return theInstance;
00096     }
00097   private:
00098     SingletonForTimings(const SingletonForTimings & other);
00099     SingletonForTimings() : accumulatedTime(0) { }
00100   };
00101 
00102 
00111   template<class Treal, class Telement = Treal>
00112     class Matrix: public MatrixHierarchicBase<Treal, Telement> {
00113     public:
00114     typedef Telement ElementType;
00115     typedef Vector<Treal, typename ElementType::VectorType> VectorType;
00116 
00117     friend class Vector<Treal, Telement>;
00118     Matrix():MatrixHierarchicBase<Treal, Telement>(){}
00119     
00120 
00121     void allocate() {
00122       assert(!this->is_empty());
00123       assert(this->is_zero());
00124       //#define MAT_USE_ALLOC_TIMER
00125 #ifdef MAT_USE_ALLOC_TIMER
00126       mat::Time theTimer; theTimer.tic();
00127 #endif
00128       this->elements = new Telement[this->nElements()];
00129 #ifdef MAT_USE_ALLOC_TIMER
00130       SingletonForTimings::instance().addTime(theTimer.toc());
00131 #endif
00132       SizesAndBlocks colSAB;
00133       SizesAndBlocks rowSAB;
00134       for (int col = 0; col < this->cols.getNBlocks(); col++) {
00135         colSAB = this->cols.getSizesAndBlocksForLowerLevel(col);
00136         for (int row = 0; row < this->rows.getNBlocks(); row++) {
00137           /* This could be improved - precompute rowSAB as well as colSAB */
00138           rowSAB = this->rows.getSizesAndBlocksForLowerLevel(row);
00139           (*this)(row,col).resetRows(rowSAB);
00140           (*this)(row,col).resetCols(colSAB);
00141         }
00142       }
00143     }
00144 
00145     /* Full matrix assigns etc */
00146     void assignFromFull(std::vector<Treal> const & fullMat);
00147     void fullMatrix(std::vector<Treal> & fullMat) const; 
00148     void syFullMatrix(std::vector<Treal> & fullMat) const;
00149     void syUpTriFullMatrix(std::vector<Treal> & fullMat) const;
00150     
00151     /* Sparse matrix assigns etc */    
00152     void assignFromSparse(std::vector<int> const & rowind, 
00153                           std::vector<int> const & colind, 
00154                           std::vector<Treal> const & values);
00155     void assignFromSparse(std::vector<int> const & rowind, 
00156                           std::vector<int> const & colind, 
00157                           std::vector<Treal> const & values,
00158                           std::vector<int> const & indexes);
00159     /* Adds values (+=) to elements */
00160     void addValues(std::vector<int> const & rowind, 
00161                    std::vector<int> const & colind, 
00162                    std::vector<Treal> const & values);
00163     void addValues(std::vector<int> const & rowind, 
00164                    std::vector<int> const & colind, 
00165                    std::vector<Treal> const & values,
00166                    std::vector<int> const & indexes);
00167 
00168     void syAssignFromSparse(std::vector<int> const & rowind, 
00169                             std::vector<int> const & colind, 
00170                             std::vector<Treal> const & values);
00171     
00172     void syAddValues(std::vector<int> const & rowind, 
00173                      std::vector<int> const & colind, 
00174                      std::vector<Treal> const & values);
00175     
00176     void getValues(std::vector<int> const & rowind, 
00177                    std::vector<int> const & colind, 
00178                    std::vector<Treal> & values) const;
00179     void getValues(std::vector<int> const & rowind, 
00180                    std::vector<int> const & colind, 
00181                    std::vector<Treal> &,
00182                    std::vector<int> const & indexes) const;
00183     void syGetValues(std::vector<int> const & rowind, 
00184                      std::vector<int> const & colind, 
00185                      std::vector<Treal> & values) const;
00186     void getAllValues(std::vector<int> & rowind, 
00187                       std::vector<int> & colind, 
00188                       std::vector<Treal> &) const;
00189     void syGetAllValues(std::vector<int> & rowind, 
00190                         std::vector<int> & colind, 
00191                         std::vector<Treal> &) const;
00192     
00193 
00194     Matrix<Treal, Telement>& 
00195     operator=(const Matrix<Treal, Telement>& mat) {
00196       MatrixHierarchicBase<Treal, Telement>::operator=(mat);
00197       return *this;
00198     } 
00199     
00200 
00201     void clear(); 
00202     ~Matrix() {
00203       clear();
00204     }
00205     void writeToFile(std::ofstream & file) const;
00206     void readFromFile(std::ifstream & file);
00207     
00208     void random();
00209     void syRandom();
00213     void randomZeroStructure(Treal probabilityBeingZero);
00214     void syRandomZeroStructure(Treal probabilityBeingZero);
00215     
00216     template<typename TRule>
00217     void setElementsByRule(TRule & rule);
00218     template<typename TRule>
00219     void sySetElementsByRule(TRule & rule);
00220     template<typename TRule>
00221     void trSetElementsByRule(TRule & rule) {
00222       // Setting elements for triangular is the same as for symmetric
00223       sySetElementsByRule(rule);
00224     }
00225     
00226     void addIdentity(Treal alpha);                  /* C += alpha * I     */
00227 
00228     static void transpose(Matrix<Treal, Telement> const & A,
00229                           Matrix<Treal, Telement> & AT);
00230 
00231     void symToNosym();
00232     void nosymToSym();
00233 
00234     /* Basic linear algebra routines */
00235 
00236     /* Set matrix to zero (k = 0) or identity (k = 1)                     */
00237     Matrix<Treal, Telement>& operator=(int const k);
00238     
00239     Matrix<Treal, Telement>& operator*=(const Treal alpha);
00240 
00241     static void gemm(const bool tA, const bool tB, const Treal alpha, 
00242                      const Matrix<Treal, Telement>& A, 
00243                      const Matrix<Treal, Telement>& B, 
00244                      const Treal beta, 
00245                      Matrix<Treal, Telement>& C);
00246     static void symm(const char side, const char uplo, const Treal alpha, 
00247                      const Matrix<Treal, Telement>& A, 
00248                      const Matrix<Treal, Telement>& B, 
00249                      const Treal beta, 
00250                      Matrix<Treal, Telement>& C);
00251     static void syrk(const char uplo, const bool tA, const Treal alpha, 
00252                      const Matrix<Treal, Telement>& A, 
00253                      const Treal beta, 
00254                      Matrix<Treal, Telement>& C);
00255     /* C = alpha * A * A + beta * C  where A and C are symmetric */
00256     static void sysq(const char uplo, const Treal alpha, 
00257                      const Matrix<Treal, Telement>& A, 
00258                      const Treal beta, 
00259                      Matrix<Treal, Telement>& C);
00260     /* C = alpha * A * B + beta * C where A and B are symmetric */
00261     static void ssmm(const Treal alpha, 
00262                      const Matrix<Treal, Telement>& A, 
00263                      const Matrix<Treal, Telement>& B, 
00264                      const Treal beta, 
00265                      Matrix<Treal, Telement>& C);
00266     /* C = alpha * A * B + beta * C where A and B are symmetric
00267      * and only the upper triangle of C is computed.
00268      */
00269     static void ssmm_upper_tr_only(const Treal alpha, 
00270                                    const Matrix<Treal, Telement>& A, 
00271                                    const Matrix<Treal, Telement>& B, 
00272                                    const Treal beta, 
00273                                    Matrix<Treal, Telement>& C);
00274 
00275     static void trmm(const char side, const char uplo, const bool tA, 
00276                      const Treal alpha, 
00277                      const Matrix<Treal, Telement>& A, 
00278                      Matrix<Treal, Telement>& B);
00279 
00280     /* Frobenius norms */
00281 
00282     /* Returns the Frobenius norm of the matrix.    */
00283     inline Treal frob() const { 
00284       return template_blas_sqrt(this->frobSquared());  
00285     } 
00286     /* Returns the squared Frobenius norm */
00287     Treal frobSquared() const;     
00288     inline Treal syFrob() const {
00289       return template_blas_sqrt(this->syFrobSquared());
00290     }
00291     Treal syFrobSquared() const;
00292     
00293     inline static Treal frobDiff
00294     (const Matrix<Treal, Telement>& A,
00295      const Matrix<Treal, Telement>& B) {
00296       return template_blas_sqrt(frobSquaredDiff(A, B));
00297     }
00298     static Treal frobSquaredDiff
00299     (const Matrix<Treal, Telement>& A,
00300      const Matrix<Treal, Telement>& B);
00301     
00302     inline static Treal syFrobDiff
00303     (const Matrix<Treal, Telement>& A,
00304      const Matrix<Treal, Telement>& B) {
00305       return template_blas_sqrt(syFrobSquaredDiff(A, B));       
00306     }
00307     static Treal syFrobSquaredDiff
00308     (const Matrix<Treal, Telement>& A,
00309      const Matrix<Treal, Telement>& B);      
00310 
00311     /* Traces */
00312     Treal trace() const;
00313     static Treal trace_ab(const Matrix<Treal, Telement>& A,
00314                           const Matrix<Treal, Telement>& B); 
00315     static Treal trace_aTb(const Matrix<Treal, Telement>& A,
00316                            const Matrix<Treal, Telement>& B); 
00317     static Treal sy_trace_ab(const Matrix<Treal, Telement>& A,
00318                              const Matrix<Treal, Telement>& B); 
00319     
00320     static void add(const Treal alpha,   /* B += alpha * A */
00321                     const Matrix<Treal, Telement>& A, 
00322                     Matrix<Treal, Telement>& B);
00323     void assign(Treal const  alpha,   /* *this = alpha * A */
00324                 Matrix<Treal, Telement> const & A);
00325 
00326 
00327     /********** Help functions for thresholding */
00328     //  int getnnzBlocksLowestLevel() const;
00329     void getFrobSqLowestLevel(std::vector<Treal> & frobsq) const;
00330     void frobThreshLowestLevel
00331     (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix);
00332 
00333     void getFrobSqElementLevel(std::vector<Treal> & frobsq) const;
00334     void frobThreshElementLevel
00335     (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix);
00336     
00337 
00338 #if 0
00339     inline void frobThreshLowestLevel
00340     (Treal const threshold, 
00341      Matrix<Treal, Telement> * ErrorMatrix) {
00342       bool a,b;
00343       frobThreshLowestLevel(threshold, ErrorMatrix, a, b);
00344     }
00345 #endif
00346 
00349     void assignFrobNormsLowestLevel
00350     ( Matrix<Treal, Matrix<Treal, Telement> > const & A ); 
00352     void syAssignFrobNormsLowestLevel
00353     ( Matrix<Treal, Matrix<Treal, Telement> > const & A ); 
00354 
00358     void assignDiffFrobNormsLowestLevel
00359     ( Matrix<Treal, Matrix<Treal, Telement> > const & A,
00360       Matrix<Treal, Matrix<Treal, Telement> > const & B ); 
00364     void syAssignDiffFrobNormsLowestLevel
00365     ( Matrix<Treal, Matrix<Treal, Telement> > const & A,
00366       Matrix<Treal, Matrix<Treal, Telement> > const & B ); 
00367 
00370     void truncateAccordingToSparsityPattern( Matrix<Treal, Matrix<Treal, Telement> > & A ) const;    
00371 
00372 
00373       /********** End of help functions for thresholding */
00374 
00375     static void gemm_upper_tr_only(const bool tA, const bool tB, 
00376                                    const Treal alpha, 
00377                                    const Matrix<Treal, Telement>& A, 
00378                                    const Matrix<Treal, Telement>& B, 
00379                                    const Treal beta, 
00380                                    Matrix<Treal, Telement>& C);
00381     static void sytr_upper_tr_only(char const side, const Treal alpha,
00382                                    Matrix<Treal, Telement>& A,
00383                                    const Matrix<Treal, Telement>& Z);
00384     static void trmm_upper_tr_only(const char side, const char uplo, 
00385                                    const bool tA, const Treal alpha, 
00386                                    const Matrix<Treal, Telement>& A, 
00387                                    Matrix<Treal, Telement>& B);
00388     static void trsytriplemm(char const side, 
00389                              const Matrix<Treal, Telement>& Z,
00390                              Matrix<Treal, Telement>& A);
00391 
00392     inline Treal frob_thresh
00393     (Treal const threshold,
00394      Matrix<Treal, Telement> * ErrorMatrix = 0) { 
00395       return template_blas_sqrt(frob_squared_thresh(threshold * threshold, 
00396                                                     ErrorMatrix)); 
00397     }                         
00403     Treal frob_squared_thresh
00404     (Treal const threshold,
00405      Matrix<Treal, Telement> * ErrorMatrix = 0); 
00411     static void syInch(const Matrix<Treal, Telement>& A,
00412                        Matrix<Treal, Telement>& Z,
00413                        const Treal threshold = 0, 
00414                        const side looking = left,
00415                        const inchversion version = unstable);
00416 
00417     void gersgorin(Treal& lmin, Treal& lmax) const; /* Computes bounds for*/
00418     /* real part of eigenvalues. The matrix must be quadratic (of course) */
00419     void sy_gersgorin(Treal& lmin, Treal& lmax) const {
00420       Matrix<Treal, Telement> tmp = (*this);
00421       tmp.symToNosym();
00422       tmp.gersgorin(lmin, lmax);
00423       return;
00424     }
00425 
00426     void add_abs_col_sums(Treal* abscolsums) const; /* Adds the absolute   */
00427     /* column sums to the abscolsums array.                                */
00428     /* abscolsums(i) += sum(abs(C(:,i))) for all i       (C == *this)      */
00429     /* Used by e.g. gersgorin eigenvalue inclusion                         */
00430     void get_diagonal(Treal* diag) const; /*Copy diagonal to the diag array*/
00431     
00432     size_t memory_usage() const; /* Returns memory usage in bytes */
00433 
00434     /* Note: we use size_t instead of int for nnz and nvalues to avoid
00435        integer overflow. */
00436     size_t nnz() const;    
00437     size_t sy_nnz() const; 
00441     inline size_t nvalues() const {
00442       return nnz();
00443     } 
00446     size_t sy_nvalues() const; 
00453     template<typename Top> 
00454     Treal syAccumulateWith(Top & op) {
00455       Treal res = 0;
00456       if (!this->is_zero()) {
00457         for (int col = 0; col < this->ncols(); col++) {
00458           for (int row = 0; row < col; row++) 
00459             res += 2 * (*this)(row, col).geAccumulateWith(op);
00460           res += (*this)(col, col).syAccumulateWith(op);
00461         }
00462       }
00463       return res;
00464     }
00466     template<typename Top>
00467     Treal geAccumulateWith(Top & op) {
00468       Treal res = 0;
00469       if (!this->is_zero()) {
00470         for (int col = 0; col < this->ncols(); col++)
00471           for (int row = 0; row < this->nrows(); row++)
00472             res += (*this)(row, col).geAccumulateWith(op);
00473       }
00474       return res;
00475     }
00476     
00477     static inline unsigned int level() {return Telement::level() + 1;}
00478 
00479     Treal maxAbsValue() const {
00480       if (this->is_zero())
00481         return 0;
00482       else {
00483         Treal maxAbsGlobal = 0;
00484         Treal maxAbsLocal  = 0;
00485         for (int ind = 0; ind < this->nElements(); ++ind) {
00486           maxAbsLocal = this->elements[ind].maxAbsValue();
00487           maxAbsGlobal = maxAbsGlobal > maxAbsLocal ? 
00488             maxAbsGlobal : maxAbsLocal; 
00489         } /* end for */
00490         return maxAbsGlobal;
00491       } 
00492     }
00493     
00494       protected:
00495       private:
00496   }; /* end class */
00497   
00498 
00499 #if 1
00500   /* Full matrix assigns etc */
00501   template<class Treal, class Telement> 
00502     void Matrix<Treal, Telement>::
00503     assignFromFull(std::vector<Treal> const & fullMat) {
00504     int nTotalRows = this->rows.getNTotalScalars();
00505     int nTotalCols = this->cols.getNTotalScalars();
00506     assert((int)fullMat.size() == nTotalRows * nTotalCols);
00507     if (this->is_zero())
00508       allocate();
00509     for (int col = 0; col < this->ncols(); col++) 
00510       for (int row = 0; row < this->nrows(); row++) 
00511         (*this)(row, col).assignFromFull(fullMat);
00512   }
00513   
00514   template<class Treal, class Telement> 
00515     void Matrix<Treal, Telement>::
00516     fullMatrix(std::vector<Treal> & fullMat) const {
00517     int nTotalRows = this->rows.getNTotalScalars();
00518     int nTotalCols = this->cols.getNTotalScalars();
00519     fullMat.resize(nTotalRows * nTotalCols);
00520     if (this->is_zero()) {
00521       int rowOffset = this->rows.getOffset();
00522       int colOffset = this->cols.getOffset();
00523       for (int col = 0; col < this->nScalarsCols(); col++)
00524         for (int row = 0; row < this->nScalarsRows(); row++) 
00525           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
00526     } 
00527     else {
00528       for (int col = 0; col < this->ncols(); col++) 
00529         for (int row = 0; row < this->nrows(); row++) 
00530           (*this)(row, col).fullMatrix(fullMat);
00531     }
00532   }
00533   
00534   template<class Treal, class Telement> 
00535     void Matrix<Treal, Telement>::
00536     syFullMatrix(std::vector<Treal> & fullMat) const {
00537     int nTotalRows = this->rows.getNTotalScalars();
00538     int nTotalCols = this->cols.getNTotalScalars();
00539     fullMat.resize(nTotalRows * nTotalCols);
00540     if (this->is_zero()) {
00541       int rowOffset = this->rows.getOffset();
00542       int colOffset = this->cols.getOffset();
00543       for (int col = 0; col < this->nScalarsCols(); col++)
00544         for (int row = 0; row <= col; row++) {
00545           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
00546           fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
00547         }
00548     } 
00549     else {
00550       for (int col = 0; col < this->ncols(); col++) {
00551         for (int row = 0; row < col; row++) 
00552           (*this)(row, col).syUpTriFullMatrix(fullMat); 
00553         (*this)(col, col).syFullMatrix(fullMat);        
00554       }
00555     }
00556   }
00557   
00558   template<class Treal, class Telement> 
00559     void Matrix<Treal, Telement>::
00560     syUpTriFullMatrix(std::vector<Treal> & fullMat) const {
00561     int nTotalRows = this->rows.getNTotalScalars();
00562     int nTotalCols = this->cols.getNTotalScalars();
00563     fullMat.resize(nTotalRows * nTotalCols);
00564     if (this->is_zero()) {
00565       int rowOffset = this->rows.getOffset();
00566       int colOffset = this->cols.getOffset();
00567       for (int col = 0; col < this->nScalarsCols(); col++)
00568         for (int row = 0; row <= this->nScalarsRows(); row++) {
00569           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
00570           fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
00571         }
00572     } 
00573     else {
00574       for (int col = 0; col < this->ncols(); col++) 
00575         for (int row = 0; row < this->nrows(); row++) 
00576           (*this)(row, col).syUpTriFullMatrix(fullMat);
00577     }
00578   }
00579   
00580 #endif
00581 
00582 
00583   template<class Treal, class Telement>
00584     void Matrix<Treal, Telement>::
00585     assignFromSparse(std::vector<int> const & rowind, 
00586                      std::vector<int> const & colind, 
00587                      std::vector<Treal> const & values) {
00588     assert(rowind.size() == colind.size() &&
00589            rowind.size() == values.size());
00590     std::vector<int> indexes(values.size());
00591     for (unsigned int ind = 0; ind < values.size(); ++ind) 
00592       indexes[ind] = ind;
00593     assignFromSparse(rowind, colind, values, indexes);
00594   }
00595   
00596   template<class Treal, class Telement>
00597     void Matrix<Treal, Telement>::
00598     assignFromSparse(std::vector<int> const & rowind, 
00599                      std::vector<int> const & colind, 
00600                      std::vector<Treal> const & values,
00601                      std::vector<int> const & indexes) {
00602     if (indexes.empty()) {
00603       this->clear();
00604       return;
00605     }
00606     if (this->is_zero())
00607       allocate();
00608     
00609     std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
00610     std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());    
00611     int currentInd = 0;
00612     
00613 
00614     std::vector<int>::const_iterator it;
00615     for ( it = indexes.begin(); it < indexes.end(); it++ )
00616       columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
00617     
00618     /* Go through all column buckets. */
00619     for (int col = 0; col < this->ncols(); col++) {
00620       /* Go through current column bucket and distribute to row buckets. */
00621       while (!columnBuckets[col].empty()) {
00622         currentInd = columnBuckets[col].back();
00623         columnBuckets[col].pop_back();
00624         rowBuckets[ this->rows.whichBlock
00625                     ( rowind[currentInd] ) ].push_back (currentInd);
00626       }
00627       /* Make calls to lower level for every row bucket. */
00628       for (int row = 0; row < this->nrows(); row++) {
00629         (*this)(row,col).assignFromSparse(rowind, colind, values, rowBuckets[row]);
00630         rowBuckets[row].clear();
00631       } /* end row loop */
00632     } /* end column loop */
00633   }  
00634 
00635   template<class Treal, class Telement>
00636     void Matrix<Treal, Telement>::
00637     addValues(std::vector<int> const & rowind, 
00638               std::vector<int> const & colind, 
00639               std::vector<Treal> const & values) {
00640     assert(rowind.size() == colind.size() &&
00641            rowind.size() == values.size());
00642     std::vector<int> indexes(values.size());
00643     for (unsigned int ind = 0; ind < values.size(); ++ind) 
00644       indexes[ind] = ind;
00645     addValues(rowind, colind, values, indexes);
00646   }
00647   
00648   template<class Treal, class Telement>
00649     void Matrix<Treal, Telement>::
00650     addValues(std::vector<int> const & rowind, 
00651               std::vector<int> const & colind, 
00652               std::vector<Treal> const & values,
00653               std::vector<int> const & indexes) {
00654     if (indexes.empty())
00655       return;
00656     if (this->is_zero())
00657       allocate();
00658     
00659     std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
00660     std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());    
00661     int currentInd = 0;
00662     
00663     std::vector<int>::const_iterator it;
00664     for ( it = indexes.begin(); it < indexes.end(); it++ )
00665       columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
00666     
00667     /* Go through all column buckets. */
00668     for (int col = 0; col < this->ncols(); col++) {
00669       /* Go through current column bucket and distribute to row buckets. */
00670       while (!columnBuckets[col].empty()) {
00671         currentInd = columnBuckets[col].back();
00672         columnBuckets[col].pop_back();
00673         rowBuckets[ this->rows.whichBlock
00674                     ( rowind[currentInd] ) ].push_back (currentInd);
00675       }
00676       /* Make calls to lower level for every row bucket. */
00677       for (int row = 0; row < this->nrows(); row++) {
00678         (*this)(row,col).addValues(rowind, colind, values, rowBuckets[row]);
00679         rowBuckets[row].clear();
00680       } /* end row loop */
00681     } /* end column loop */
00682   }
00683 
00684   template<class Treal, class Telement>
00685     void Matrix<Treal, Telement>::
00686     syAssignFromSparse(std::vector<int> const & rowind, 
00687                        std::vector<int> const & colind, 
00688                        std::vector<Treal> const & values) {
00689     assert(rowind.size() == colind.size() &&
00690            rowind.size() == values.size());
00691     bool upperTriangleOnly = true;
00692     for (unsigned int ind = 0; ind < values.size(); ++ind) {
00693       upperTriangleOnly = 
00694         upperTriangleOnly && colind[ind] >= rowind[ind];
00695     }
00696     if (!upperTriangleOnly)
00697       throw Failure("Matrix<Treal, Telement>::"
00698                     "syAddValues(std::vector<int> const &, "
00699                     "std::vector<int> const &, "
00700                     "std::vector<Treal> const &, int const): "
00701                     "Only upper triangle can contain elements when assigning "
00702                     "symmetric or triangular matrix ");
00703     assignFromSparse(rowind, colind, values);
00704   }
00705   
00706   template<class Treal, class Telement>
00707     void Matrix<Treal, Telement>::
00708     syAddValues(std::vector<int> const & rowind, 
00709                 std::vector<int> const & colind, 
00710                 std::vector<Treal> const & values) {
00711     assert(rowind.size() == colind.size() &&
00712            rowind.size() == values.size());
00713     bool upperTriangleOnly = true;
00714     for (unsigned int ind = 0; ind < values.size(); ++ind) {
00715       upperTriangleOnly = 
00716         upperTriangleOnly && colind[ind] >= rowind[ind];
00717     }
00718     if (!upperTriangleOnly)
00719       throw Failure("Matrix<Treal, Telement>::"
00720                     "syAddValues(std::vector<int> const &, "
00721                     "std::vector<int> const &, "
00722                     "std::vector<Treal> const &, int const): "
00723                     "Only upper triangle can contain elements when assigning "
00724                     "symmetric or triangular matrix ");
00725     addValues(rowind, colind, values);
00726   }
00727 
00728   template<class Treal, class Telement>
00729     void Matrix<Treal, Telement>::
00730     getValues(std::vector<int> const & rowind, 
00731               std::vector<int> const & colind, 
00732               std::vector<Treal> & values) const {
00733     assert(rowind.size() == colind.size());
00734     values.resize(rowind.size(),0);
00735     std::vector<int> indexes(rowind.size());
00736     for (unsigned int ind = 0; ind < rowind.size(); ++ind) 
00737       indexes[ind] = ind;
00738     getValues(rowind, colind, values, indexes);
00739   }
00740   
00741   template<class Treal, class Telement>
00742     void Matrix<Treal, Telement>::
00743     getValues(std::vector<int> const & rowind, 
00744               std::vector<int> const & colind, 
00745               std::vector<Treal> & values,
00746               std::vector<int> const & indexes) const {
00747     assert(!this->is_empty());
00748     if (indexes.empty())
00749       return;
00750     std::vector<int>::const_iterator it;
00751     if (this->is_zero()) {
00752       for ( it = indexes.begin(); it < indexes.end(); it++ )
00753         values[*it] = 0;
00754       return;
00755     }
00756     
00757     std::vector<std::vector<int> > columnBuckets (this->cols.getNBlocks());
00758     std::vector<std::vector<int> > rowBuckets (this->rows.getNBlocks());    
00759     int currentInd = 0;
00760     
00761     for ( it = indexes.begin(); it < indexes.end(); it++ )
00762       columnBuckets[ this->cols.whichBlock( colind[*it] ) ].push_back (*it);
00763     
00764     /* Go through all column buckets. */
00765     for (int col = 0; col < this->ncols(); col++) {
00766       /* Go through current column bucket and distribute to row buckets. */
00767       while (!columnBuckets[col].empty()) {
00768         currentInd = columnBuckets[col].back();
00769         columnBuckets[col].pop_back();
00770         rowBuckets[ this->rows.whichBlock
00771                     ( rowind[currentInd] ) ].push_back (currentInd);
00772       }
00773       /* Make calls to lower level for every row bucket. */
00774       for (int row = 0; row < this->nrows(); row++) {
00775         (*this)(row,col).getValues(rowind, colind, values, rowBuckets[row]);
00776         rowBuckets[row].clear();
00777       } /* end row loop */
00778     } /* end column loop */
00779   }
00780   
00781   template<class Treal, class Telement>
00782     void Matrix<Treal, Telement>::
00783     syGetValues(std::vector<int> const & rowind, 
00784                 std::vector<int> const & colind, 
00785                 std::vector<Treal> & values) const {
00786     assert(rowind.size() == colind.size());
00787     bool upperTriangleOnly = true;
00788     for (int unsigned ind = 0; ind < rowind.size(); ++ind) {
00789       upperTriangleOnly = 
00790         upperTriangleOnly && colind[ind] >= rowind[ind];
00791     }
00792     if (!upperTriangleOnly)
00793       throw Failure("Matrix<Treal, Telement>::"
00794                     "syGetValues(std::vector<int> const &, "
00795                     "std::vector<int> const &, "
00796                     "std::vector<Treal> const &, int const): "
00797                     "Only upper triangle when retrieving elements from "
00798                     "symmetric or triangular matrix ");
00799     getValues(rowind, colind, values);
00800   }
00801   
00802 
00803   template<class Treal, class Telement>
00804     void Matrix<Treal, Telement>::
00805     getAllValues(std::vector<int> & rowind, 
00806                  std::vector<int> & colind, 
00807                  std::vector<Treal> & values) const {
00808     assert(rowind.size() == colind.size() &&
00809            rowind.size() == values.size());
00810     if (!this->is_zero()) 
00811       for (int ind = 0; ind < this->nElements(); ++ind)
00812         this->elements[ind].getAllValues(rowind, colind, values);
00813   }
00814 
00815   template<class Treal, class Telement>
00816     void Matrix<Treal, Telement>::
00817     syGetAllValues(std::vector<int> & rowind, 
00818                    std::vector<int> & colind, 
00819                    std::vector<Treal> & values) const {
00820     assert(rowind.size() == colind.size() &&
00821            rowind.size() == values.size());
00822     if (!this->is_zero()) 
00823       for (int col = 0; col < this->ncols(); ++col) {
00824         for (int row = 0; row < col; ++row)
00825           (*this)(row, col).getAllValues(rowind, colind, values);
00826         (*this)(col, col).syGetAllValues(rowind, colind, values);
00827       }
00828   }
00829   
00830   
00831   template<class Treal, class Telement>
00832     void Matrix<Treal, Telement>::clear() {
00833     if (!this->is_zero())
00834       for (int i = 0; i < this->nElements(); i++) 
00835         this->elements[i].clear();
00836     delete[] this->elements;
00837     this->elements = 0;
00838   }
00839   
00840   template<class Treal, class Telement>
00841     void Matrix<Treal, Telement>:: 
00842     writeToFile(std::ofstream & file) const {
00843     int const ZERO = 0;
00844     int const ONE  = 1;
00845     if (this->is_zero()) {
00846       char * tmp = (char*)&ZERO;
00847       file.write(tmp,sizeof(int));
00848     }
00849     else {
00850       char * tmp = (char*)&ONE;
00851       file.write(tmp,sizeof(int));
00852       for (int i = 0; i < this->nElements(); i++)
00853         this->elements[i].writeToFile(file);
00854     }
00855   }
00856   template<class Treal, class Telement>
00857     void Matrix<Treal, Telement>:: 
00858     readFromFile(std::ifstream & file) {
00859     int const ZERO = 0;
00860     int const ONE  = 1;
00861     char tmp[sizeof(int)];
00862     file.read(tmp, (std::ifstream::pos_type)sizeof(int));
00863     switch ((int)*tmp) {
00864     case ZERO:
00865       this->clear();
00866       break;
00867     case ONE:
00868       if (this->is_zero()) 
00869         allocate();
00870       for (int i = 0; i < this->nElements(); i++)
00871         this->elements[i].readFromFile(file);
00872       break;
00873     default:
00874       throw Failure("Matrix<Treal, Telement>::" 
00875                     "readFromFile(std::ifstream & file):"
00876                     "File corruption int value not 0 or 1");
00877     }
00878   }
00879 
00880   template<class Treal, class Telement> 
00881     void Matrix<Treal, Telement>::random() {
00882     if (this->is_zero()) 
00883       allocate();
00884     for (int ind = 0; ind < this->nElements(); ind++)
00885       this->elements[ind].random();     
00886   }
00887 
00888   template<class Treal, class Telement> 
00889     void Matrix<Treal, Telement>::syRandom() {
00890     if (this->is_zero()) 
00891       allocate();
00892     /* Above diagonal */
00893     for (int col = 1; col < this->ncols(); col++)
00894       for (int row = 0; row < col; row++)
00895         (*this)(row, col).random();
00896     /* Diagonal */
00897     for (int rc = 0; rc < this->nrows(); rc++)
00898       (*this)(rc,rc).syRandom();
00899   }
00900 
00901   template<class Treal, class Telement> 
00902     void Matrix<Treal, Telement>::
00903     randomZeroStructure(Treal probabilityBeingZero) {
00904     if (!this->highestLevel() && 
00905         probabilityBeingZero > rand() / (Treal)RAND_MAX)
00906       this->clear();
00907     else {
00908       if (this->is_zero()) 
00909         allocate();
00910       for (int ind = 0; ind < this->nElements(); ind++)
00911         this->elements[ind].randomZeroStructure(probabilityBeingZero);          
00912     }      
00913   }
00914 
00915   template<class Treal, class Telement> 
00916     void  Matrix<Treal, Telement>::
00917     syRandomZeroStructure(Treal probabilityBeingZero) {
00918     if (!this->highestLevel() && 
00919         probabilityBeingZero > rand() / (Treal)RAND_MAX)
00920       this->clear();
00921     else {
00922       if (this->is_zero()) 
00923         allocate();
00924       /* Above diagonal */
00925       for (int col = 1; col < this->ncols(); col++)
00926         for (int row = 0; row < col; row++)
00927           (*this)(row, col).randomZeroStructure(probabilityBeingZero);
00928       /* Diagonal */
00929       for (int rc = 0; rc < this->nrows(); rc++)
00930         (*this)(rc,rc).syRandomZeroStructure(probabilityBeingZero);
00931     }
00932   }
00933 
00934   template<class Treal, class Telement> 
00935     template<typename TRule>
00936     void Matrix<Treal, Telement>::
00937     setElementsByRule(TRule & rule) {
00938     if (this->is_zero()) 
00939       allocate();
00940     for (int ind = 0; ind < this->nElements(); ind++)
00941       this->elements[ind].setElementsByRule(rule);      
00942   }
00943   template<class Treal, class Telement> 
00944     template<typename TRule>
00945     void Matrix<Treal, Telement>::
00946     sySetElementsByRule(TRule & rule) {
00947     if (this->is_zero()) 
00948       allocate();
00949     /* Above diagonal */
00950     for (int col = 1; col < this->ncols(); col++)
00951       for (int row = 0; row < col; row++)
00952         (*this)(row, col).setElementsByRule(rule);
00953     /* Diagonal */
00954     for (int rc = 0; rc < this->nrows(); rc++)
00955       (*this)(rc,rc).sySetElementsByRule(rule);
00956   }
00957   
00958 
00959   template<class Treal, class Telement> 
00960     void Matrix<Treal, Telement>::
00961     addIdentity(Treal alpha) {
00962     if (this->is_empty())
00963       throw Failure("Matrix<Treal, Telement>::addIdentity(Treal): "
00964                     "Cannot add identity to empty matrix.");
00965     if (this->ncols() != this->nrows()) 
00966       throw Failure("Matrix<Treal, Telement>::addIdentity(Treal): "
00967                     "Matrix must be square to add identity");
00968     if (this->is_zero()) 
00969       allocate();
00970     for (int ind = 0; ind < this->ncols(); ind++)
00971       (*this)(ind,ind).addIdentity(alpha);
00972   }
00973   
00974   template<class Treal, class Telement> 
00975     void Matrix<Treal, Telement>::
00976     transpose(Matrix<Treal, Telement> const & A,
00977               Matrix<Treal, Telement> & AT) {
00978     if (A.is_zero()) { /* Condition also matches empty matrices. */
00979       AT.rows = A.cols;
00980       AT.cols = A.rows;
00981       delete[] AT.elements;
00982       AT.elements = 0;
00983       return;
00984     }
00985     if (AT.is_zero() || (AT.nElements() != A.nElements())) {
00986       delete[] AT.elements;
00987       AT.elements = new Telement[A.nElements()];
00988     }
00989     AT.cols = A.rows;
00990     AT.rows = A.cols;
00991     for (int row = 0; row < AT.nrows(); row++)
00992       for (int col = 0; col < AT.ncols(); col++)
00993         Telement::transpose(A(col,row), AT(row,col));
00994   }
00995 
00996   template<class Treal, class Telement> 
00997     void Matrix<Treal, Telement>::
00998     symToNosym() {
00999     try {
01000       if (this->nrows() == this->ncols()) {
01001         if (!this->is_zero()) {
01002           /* Fix the diagonal: */
01003           for (int rc = 0; rc < this->ncols(); rc++)
01004             (*this)(rc, rc).symToNosym();
01005           /* Fix the lower triangle */
01006           for (int row = 1; row < this->nrows(); row++)
01007             for (int col = 0; col < row; col++)
01008               Telement::transpose((*this)(col, row), (*this)(row, col));
01009         }
01010       }
01011       else
01012         throw Failure("Matrix<Treal, Telement>::symToNosym(): "
01013                       "Only quadratic matrices can be symmetric");      
01014     }
01015     catch(Failure& f) {
01016       std::cout<<"Failure caught:Matrix<Treal, Telement>::symToNosym()"
01017                <<std::endl;
01018       throw f;
01019     }
01020   }
01021 
01022   template<class Treal, class Telement> 
01023     void Matrix<Treal, Telement>::
01024     nosymToSym() {
01025     if (this->nrows() == this->ncols()) {
01026       if (!this->is_zero()) {
01027         /* Fix the diagonal: */
01028         for (int rc = 0; rc < this->ncols(); rc++)
01029           (*this)(rc, rc).nosymToSym();
01030         /* Fix the lower triangle */
01031         for (int col = 0; col < this->ncols() - 1; col++)
01032           for (int row = col + 1; row < this->nrows(); row++)
01033             (*this)(row, col).clear();
01034       }    
01035     }
01036     else
01037       throw Failure("Matrix<Treal, Telement>::nosymToSym(): "
01038                     "Only quadratic matrices can be symmetric");      
01039   }
01040   
01041   /* Basic linear algebra routines */
01042 
01043   template<class Treal, class Telement>
01044     Matrix<Treal, Telement>& 
01045     Matrix<Treal, Telement>::operator=(int const k) {
01046     switch (k) {
01047     case 0:
01048       this->clear();
01049       break;
01050     case 1:
01051       if (this->ncols() != this->nrows()) 
01052         throw Failure
01053           ("Matrix::operator=(int k = 1): "
01054            "Matrix must be quadratic to become identity matrix.");
01055       this->clear();
01056       this->allocate();
01057       for (int rc = 0; rc < this->ncols(); rc++) /*Set diagonal to identity*/
01058         (*this)(rc,rc) = 1;
01059       break;
01060     default:
01061       throw Failure("Matrix::operator=(int k) only "
01062                     "implemented for k = 0, k = 1");
01063     }
01064     return *this;
01065   }
01066 
01067 
01068   template<class Treal, class Telement>
01069     Matrix<Treal, Telement>& Matrix<Treal, Telement>:: 
01070     operator*=(const Treal alpha) {
01071     if (!this->is_zero() && alpha != 1) {
01072       for (int ind = 0; ind < this->nElements(); ind++)
01073         this->elements[ind] *= alpha;
01074     }
01075     return *this;
01076   }
01077 
01078   /* C = beta * C + alpha * A * B */
01079   template<class Treal, class Telement>
01080     void Matrix<Treal, Telement>::
01081     gemm(const bool tA, const bool tB, const Treal alpha, 
01082          const Matrix<Treal, Telement>& A, 
01083          const Matrix<Treal, Telement>& B, 
01084          const Treal beta, 
01085          Matrix<Treal, Telement>& C) {
01086     assert(!A.is_empty());
01087     assert(!B.is_empty());
01088     if (C.is_empty()) {
01089       assert(beta == 0);
01090       if (tA)
01091         C.resetRows(A.cols);
01092       else
01093         C.resetRows(A.rows);
01094       if (tB)
01095         C.resetCols(B.rows);
01096       else
01097         C.resetCols(B.cols);
01098     }      
01099     
01100     if ((A.is_zero() || B.is_zero() || alpha == 0) && 
01101         (C.is_zero() || beta == 0))
01102       C = 0;
01103     else {
01104       Treal beta_tmp = beta;
01105       if (C.is_zero()) {
01106         C.allocate();
01107         beta_tmp = 0;
01108       }
01109       if (!A.is_zero() && !B.is_zero() && alpha != 0) {
01110         MAT_OMP_INIT;
01111         if (!tA && !tB) 
01112           if (A.ncols() == B.nrows() && 
01113               A.nrows() == C.nrows() &&
01114               B.ncols() == C.ncols()) { 
01115 #ifdef _OPENMP
01116 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 
01117 #endif
01118             for (int col = 0; col < C.ncols(); col++) {
01119               MAT_OMP_START;
01120               for (int row = 0; row < C.nrows(); row++) {
01121                 Telement::gemm(tA, tB, alpha, 
01122                                A(row, 0), B(0, col), 
01123                                beta_tmp, 
01124                                C(row, col));
01125                 for(int cola = 1; cola < A.ncols(); cola++)
01126                   Telement::gemm(tA, tB, alpha, 
01127                                  A(row, cola), B(cola, col), 
01128                                  1.0, 
01129                                  C(row, col));
01130               }
01131               MAT_OMP_END;
01132             } /* end omp for */
01133           }
01134           else 
01135             throw Failure("Matrix<class Treal, class Telement>::"
01136                           "gemm(bool, bool, Treal, "
01137                           "Matrix<Treal, Telement>, "
01138                           "Matrix<Treal, Telement>, Treal, "
01139                           "Matrix<Treal, Telement>): "
01140                           "Incorrect matrixdimensions for multiplication");
01141         else if (tA && !tB)
01142           if (A.nrows() == B.nrows() &&
01143               A.ncols() == C.nrows() &&
01144               B.ncols() == C.ncols()) {
01145 #ifdef _OPENMP
01146 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 
01147 #endif
01148             for (int col = 0; col < C.ncols(); col++) {
01149               MAT_OMP_START;
01150               for (int row = 0; row < C.nrows(); row++) {
01151                 Telement::gemm(tA, tB, alpha,
01152                                A(0,row), B(0,col),
01153                                beta_tmp,
01154                                C(row,col));
01155                 for(int cola = 1; cola < A.nrows(); cola++)
01156                   Telement::gemm(tA, tB, alpha, 
01157                                  A(cola, row), B(cola, col), 
01158                                  1.0, 
01159                                  C(row,col));
01160               }
01161               MAT_OMP_END;
01162             } /* end omp for */
01163           }
01164           else
01165             throw Failure("Matrix<class Treal, class Telement>::"
01166                           "gemm(bool, bool, Treal, "
01167                           "Matrix<Treal, Telement>, "
01168                           "Matrix<Treal, Telement>, Treal, "
01169                           "Matrix<Treal, Telement>): "
01170                           "Incorrect matrixdimensions for multiplication");
01171         else if (!tA && tB)
01172           if (A.ncols() == B.ncols() && 
01173               A.nrows() == C.nrows() &&
01174               B.nrows() == C.ncols()) {
01175 #ifdef _OPENMP
01176 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 
01177 #endif
01178             for (int col = 0; col < C.ncols(); col++) {
01179               MAT_OMP_START;
01180               for (int row = 0; row < C.nrows(); row++) {
01181                 Telement::gemm(tA, tB, alpha, 
01182                                A(row, 0), B(col, 0), 
01183                                beta_tmp, 
01184                                C(row,col));
01185                 for(int cola = 1; cola < A.ncols(); cola++)
01186                   Telement::gemm(tA, tB, alpha, 
01187                                  A(row, cola), B(col, cola), 
01188                                  1.0, 
01189                                  C(row,col));
01190               }
01191               MAT_OMP_END;
01192             } /* end omp for */
01193           }
01194           else 
01195             throw Failure("Matrix<class Treal, class Telement>::"
01196                           "gemm(bool, bool, Treal, "
01197                           "Matrix<Treal, Telement>, "
01198                           "Matrix<Treal, Telement>, Treal, "
01199                           "Matrix<Treal, Telement>): "
01200                           "Incorrect matrixdimensions for multiplication");
01201         else if (tA && tB)
01202           if (A.nrows() == B.ncols() && 
01203               A.ncols() == C.nrows() &&
01204               B.nrows() == C.ncols()) {
01205 #ifdef _OPENMP
01206 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 
01207 #endif
01208             for (int col = 0; col < C.ncols(); col++) {
01209               MAT_OMP_START;
01210               for (int row = 0; row < C.nrows(); row++) {
01211                 Telement::gemm(tA, tB, alpha, 
01212                                A(0, row), B(col, 0), 
01213                                beta_tmp, 
01214                                C(row,col));
01215                 for(int cola = 1; cola < A.nrows(); cola++)
01216                   Telement::gemm(tA, tB, alpha, 
01217                                  A(cola, row), B(col, cola), 
01218                                  1.0, 
01219                                  C(row,col));
01220               }
01221               MAT_OMP_END;
01222             } /* end omp for */
01223           }
01224           else 
01225             throw Failure("Matrix<class Treal, class Telement>::"
01226                           "gemm(bool, bool, Treal, "
01227                           "Matrix<Treal, Telement>, "
01228                           "Matrix<Treal, Telement>, "
01229                           "Treal, Matrix<Treal, Telement>): "
01230                           "Incorrect matrixdimensions for multiplication");
01231         else throw Failure("Matrix<class Treal, class Telement>::"
01232                            "gemm(bool, bool, Treal, "
01233                            "Matrix<Treal, Telement>, "
01234                            "Matrix<Treal, Telement>, Treal, "
01235                            "Matrix<Treal, Telement>):"
01236                            "Very strange error!!");
01237         MAT_OMP_FINALIZE;
01238       }    
01239       else
01240         C *= beta;
01241     }
01242   }
01243 
01244   template<class Treal, class Telement>
01245     void Matrix<Treal, Telement>::
01246     symm(const char side, const char uplo, const Treal alpha, 
01247          const Matrix<Treal, Telement>& A, 
01248          const Matrix<Treal, Telement>& B, 
01249          const Treal beta, 
01250          Matrix<Treal, Telement>& C) {
01251     assert(!A.is_empty());
01252     assert(!B.is_empty());
01253     assert(A.nrows() == A.ncols());
01254     int dimA = A.nrows();
01255     if (C.is_empty()) {
01256       assert(beta == 0);
01257       if (side =='L') {
01258         C.resetRows(A.rows);
01259         C.resetCols(B.cols);
01260       } 
01261       else {
01262         assert(side == 'R');
01263         C.resetRows(B.rows);
01264         C.resetCols(A.cols);
01265       }
01266     }      
01267     if ((A.is_zero() || B.is_zero() || alpha == 0) && 
01268         (C.is_zero() || beta == 0))
01269       C = 0;
01270     else {
01271       if (uplo == 'U') {
01272         Treal beta_tmp = beta;
01273         if (C.is_zero()) {
01274           C.allocate();
01275           beta_tmp = 0;
01276         }
01277         if (!A.is_zero() && !B.is_zero() && alpha != 0) {
01278           MAT_OMP_INIT;
01279           if (side =='L') 
01280             if (dimA == B.nrows() && 
01281                 dimA == C.nrows() &&
01282                 B.ncols() == C.ncols()) {
01283 #ifdef _OPENMP
01284 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 
01285 #endif
01286               for (int col = 0; col < C.ncols(); col++) {
01287                 MAT_OMP_START;
01288                 for (int row = 0; row < C.nrows(); row++) {
01289                   /* Diagonal element in matrix A */
01290                   Telement::symm(side, uplo, alpha, A(row, row),
01291                                  B(row, col), beta_tmp, C(row, col));
01292                   /* Elements below diagonal in A */
01293                   for(int ind = 0; ind < row; ind++)
01294                     Telement::gemm(true, false, alpha, 
01295                                    A(ind, row), B(ind, col), 
01296                                    1.0, C(row,col));
01297                   /* Elements above diagonal in A */
01298                   for(int ind = row + 1; ind < dimA; ind++)
01299                     Telement::gemm(false, false, alpha, 
01300                                    A(row, ind), B(ind, col), 
01301                                    1.0, C(row,col));
01302                 }
01303                 MAT_OMP_END;
01304               } /* end omp for */
01305             }
01306             else 
01307               throw Failure("Matrix<class Treal, class Telement>"
01308                             "::symm(bool, bool, Treal, Matrix<Treal, "
01309                             "Telement>, Matrix<Treal, Telement>, "
01310                             "Treal, Matrix<Treal, Telement>): "
01311                             "Incorrect matrixdimensions for multiplication");
01312           else { /* side == 'R' */
01313             assert(side == 'R');
01314             if (B.ncols() == dimA && 
01315                 B.nrows() == C.nrows() &&
01316                 dimA == C.ncols()) {
01317 #ifdef _OPENMP
01318 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 
01319 #endif
01320               for (int col = 0; col < C.ncols(); col++) {
01321                 MAT_OMP_START;
01322                 for (int row = 0; row < C.nrows(); row++) {
01323                   /* Diagonal element in matrix A */
01324                   Telement::symm(side, uplo, alpha, A(col, col),
01325                                  B(row, col), beta_tmp, C(row, col));
01326                   /* Elements below diagonal in A */
01327                   for(int ind = col + 1; ind < dimA; ind++)
01328                     Telement::gemm(false, true, alpha, 
01329                                    B(row, ind), A(col, ind),  
01330                                    1.0, C(row,col));
01331                   /* Elements above diagonal in A */
01332                   for(int ind = 0; ind < col; ind++)
01333                     Telement::gemm(false, false, alpha, 
01334                                    B(row, ind), A(ind, col),  
01335                                    1.0, C(row,col));
01336                 }
01337                 MAT_OMP_END;
01338               } /* end omp for */
01339             }
01340             else 
01341               throw Failure("Matrix<class Treal, class Telement>"
01342                             "::symm(bool, bool, Treal, Matrix<Treal, "
01343                             "Telement>, Matrix<Treal, Telement>, "
01344                             "Treal, Matrix<Treal, Telement>): "
01345                             "Incorrect matrixdimensions for multiplication");
01346           }
01347           MAT_OMP_FINALIZE;
01348         }
01349         else
01350           C *= beta;
01351       }
01352       else
01353         throw Failure("Matrix<class Treal, class Telement>::"
01354                       "symm only implemented for symmetric matrices in "
01355                       "upper triangular storage");
01356     }    
01357   }
01358 
01359   template<class Treal, class Telement>
01360     void Matrix<Treal, Telement>:: 
01361     syrk(const char uplo, const bool tA, const Treal alpha, 
01362          const Matrix<Treal, Telement>& A, 
01363          const Treal beta, 
01364          Matrix<Treal, Telement>& C) {
01365     assert(!A.is_empty());
01366     if (C.is_empty()) {
01367       assert(beta == 0);
01368       if (tA) {
01369         C.resetRows(A.cols);
01370         C.resetCols(A.cols);
01371       } 
01372       else {
01373         C.resetRows(A.rows);
01374         C.resetCols(A.rows);
01375       }
01376     }
01377     
01378     if (C.nrows() == C.ncols() && 
01379         ((!tA && A.nrows() == C.nrows()) || (tA && A.ncols() == C.nrows())))
01380       if (alpha != 0 && !A.is_zero()) {
01381         Treal beta_tmp = beta;
01382         if (C.is_zero()) {
01383           C.allocate();
01384           beta_tmp = 0;
01385         }
01386         MAT_OMP_INIT;
01387         if (!tA && uplo == 'U') { /* C = alpha * A * A' + beta * C */
01388 #ifdef _OPENMP
01389 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
01390 #endif
01391           {
01392 #ifdef _OPENMP
01393 #pragma omp for  schedule(dynamic) nowait
01394 #endif
01395             for (int rc = 0; rc < C.ncols(); rc++) {
01396               MAT_OMP_START;
01397               Telement::syrk(uplo, tA, alpha, A(rc, 0), beta_tmp, C(rc, rc));
01398               for (int cola = 1; cola < A.ncols(); cola++)
01399                 Telement::syrk(uplo, tA, alpha, A(rc, cola), 1.0, C(rc, rc));
01400               MAT_OMP_END;
01401             }
01402 #ifdef _OPENMP
01403 #pragma omp for  schedule(dynamic) nowait
01404 #endif
01405             for (int row = 0; row < C.nrows() - 1; row++) {
01406               MAT_OMP_START;
01407               for (int col = row + 1; col < C.ncols(); col++) {
01408                 Telement::gemm(tA, !tA, alpha, 
01409                                A(row, 0), A(col,0), beta_tmp, C(row,col));
01410                 for (int ind = 1; ind < A.ncols(); ind++)
01411                   Telement::gemm(tA, !tA, alpha, 
01412                                  A(row, ind), A(col,ind), 1.0, C(row,col));
01413               }
01414               MAT_OMP_END;
01415             }
01416           } /* end omp parallel */
01417         } /* end if (!tA) */
01418         else if (tA && uplo == 'U') { /* C = alpha * A' * A + beta * C */
01419 #ifdef _OPENMP
01420 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
01421 #endif
01422           {
01423 #ifdef _OPENMP
01424 #pragma omp for  schedule(dynamic) nowait
01425 #endif
01426             for (int rc = 0; rc < C.ncols(); rc++) {
01427               MAT_OMP_START;
01428               Telement::syrk(uplo, tA, alpha, A(0, rc), beta_tmp, C(rc, rc));
01429               for (int rowa = 1; rowa < A.nrows(); rowa++)
01430                 Telement::syrk(uplo, tA, alpha, A(rowa, rc), 1.0, C(rc, rc));
01431               MAT_OMP_END;
01432             }
01433 #ifdef _OPENMP
01434 #pragma omp for  schedule(dynamic) nowait
01435 #endif
01436             for (int row = 0; row < C.nrows() - 1; row++) {
01437               MAT_OMP_START;
01438               for (int col = row + 1; col < C.ncols(); col++) {
01439                 Telement::gemm(tA, !tA, alpha, 
01440                                A(0, row), A(0, col), beta_tmp, C(row,col));
01441                 for (int ind = 1; ind < A.nrows(); ind++)
01442                   Telement::gemm(tA, !tA, alpha, 
01443                                  A(ind, row), A(ind, col), 1.0, C(row,col));
01444               }
01445               MAT_OMP_END;
01446             }
01447           } /* end omp parallel */
01448         } /* end if (tA) */  
01449         else
01450           throw Failure("Matrix<class Treal, class Telement>::"
01451                         "syrk not implemented for lower triangular storage");
01452         MAT_OMP_FINALIZE;
01453       }
01454       else {
01455         C *= beta;
01456       }
01457     else
01458       throw Failure("Matrix<class Treal, class Telement>::syrk: "
01459                     "Incorrect matrix dimensions for symmetric rank-k update");
01460   }
01461 
01462   template<class Treal, class Telement>
01463     void Matrix<Treal, Telement>:: 
01464     sysq(const char uplo, const Treal alpha, 
01465          const Matrix<Treal, Telement>& A, 
01466          const Treal beta, 
01467          Matrix<Treal, Telement>& C) {
01468     assert(!A.is_empty());
01469     if (C.is_empty()) {
01470       assert(beta == 0);
01471       C.resetRows(A.rows);
01472       C.resetCols(A.cols);
01473     }
01474     if (C.nrows() == C.ncols() && 
01475         A.nrows() == C.nrows() && A.nrows() == A.ncols())
01476       if (alpha != 0 && !A.is_zero()) {
01477         if (uplo == 'U') {
01478           Treal beta_tmp = beta;
01479           if (C.is_zero()) {
01480             C.allocate();
01481             beta_tmp = 0;
01482           }
01483           MAT_OMP_INIT;
01484 #ifdef _OPENMP
01485 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
01486 #endif
01487           {
01488 #ifdef _OPENMP
01489 #pragma omp for  schedule(dynamic) nowait
01490 #endif
01491             for (int rc = 0; rc < C.ncols(); rc++) {
01492               MAT_OMP_START;
01493               Telement::sysq(uplo, alpha, A(rc, rc), beta_tmp, C(rc, rc));
01494               for (int cola = 0; cola < rc; cola++)
01495                 Telement::syrk(uplo, true, alpha, A(cola, rc), 1.0, C(rc, rc));
01496               for (int cola = rc + 1; cola < A.ncols(); cola++)
01497                 Telement::syrk(uplo, false, alpha, A(rc, cola), 1.0, C(rc, rc));
01498               MAT_OMP_END;
01499             }
01500             /* Maste anvanda symm? */
01501 #ifdef _OPENMP
01502 #pragma omp for schedule(dynamic) nowait
01503 #endif
01504             for (int row = 0; row < C.nrows() - 1; row++) {
01505               MAT_OMP_START;
01506               for (int col = row + 1; col < C.ncols(); col++) {
01507                 /* First the two operations involving diagonal elements */
01508                 Telement::symm('L', 'U', alpha, A(row, row), A(row, col), 
01509                                beta_tmp, C(row, col));
01510                 Telement::symm('R', 'U', alpha, A(col, col), A(row, col), 
01511                                1.0, C(row, col));
01512                 /* First element in product is below  the diagonal */
01513                 for (int ind = 0; ind < row; ind++)
01514                   Telement::gemm(true, false, alpha, 
01515                                  A(ind, row), A(ind, col), 1.0, C(row, col));
01516                 /* None of the elements are below the diagonal     */
01517                 for (int ind = row + 1; ind < col; ind++)
01518                   Telement::gemm(false, false, alpha, 
01519                                  A(row, ind), A(ind, col), 1.0, C(row, col));
01520                 /* Second element is below the diagonal            */
01521                 for (int ind = col + 1; ind < A.ncols(); ind++)
01522                   Telement::gemm(false, true, alpha, 
01523                                  A(row, ind), A(col, ind), 1.0, C(row, col));
01524               }
01525               MAT_OMP_END;
01526             } /* end omp for */
01527           } /* end omp parallel */
01528           MAT_OMP_FINALIZE;
01529         }
01530         else
01531           throw Failure("Matrix<class Treal, class Telement>::"
01532                         "sysq only implemented for symmetric matrices in "
01533                         "upper triangular storage");;
01534       }
01535       else {
01536         C *= beta;
01537       }
01538     else
01539       throw Failure("Matrix<class Treal, class Telement>::sysq: "
01540                     "Incorrect matrix dimensions for symmetric square "
01541                     "operation");
01542   }
01543 
01544   /* C = alpha * A * B + beta * C where A and B are symmetric */
01545   template<class Treal, class Telement>
01546     void Matrix<Treal, Telement>:: 
01547     ssmm(const Treal alpha, 
01548          const Matrix<Treal, Telement>& A, 
01549          const Matrix<Treal, Telement>& B, 
01550          const Treal beta, 
01551          Matrix<Treal, Telement>& C) {
01552     assert(!A.is_empty());
01553     assert(!B.is_empty());
01554     if (C.is_empty()) {
01555       assert(beta == 0);
01556       C.resetRows(A.rows);
01557       C.resetCols(B.cols);
01558     }
01559     if (A.ncols() != B.nrows() || 
01560         A.nrows() != C.nrows() ||
01561         B.ncols() != C.ncols() ||
01562         A.nrows() != A.ncols() ||
01563         B.nrows() != B.ncols()) {
01564       throw Failure("Matrix<class Treal, class Telement>::"
01565                     "ssmm(Treal, "
01566                     "Matrix<Treal, Telement>, "
01567                     "Matrix<Treal, Telement>, Treal, "
01568                     "Matrix<Treal, Telement>): "
01569                     "Incorrect matrixdimensions for ssmm multiplication");
01570     }
01571     if ((A.is_zero() || B.is_zero() || alpha == 0) && 
01572         (C.is_zero() || beta == 0)) {
01573       C = 0;
01574       return;
01575     }
01576     Treal beta_tmp = beta;
01577     if (C.is_zero()) {
01578       C.allocate();
01579       beta_tmp = 0;
01580     }
01581     if (A.is_zero() || B.is_zero() || alpha == 0) {
01582       C *= beta;
01583       return;
01584     }
01585     MAT_OMP_INIT;
01586 #ifdef _OPENMP
01587 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 
01588 #endif
01589     for (int col = 0; col < C.ncols(); col++) {
01590       MAT_OMP_START;
01591       /* Diagonal */
01592       /* C(col, col) = alpha * A(col, col) * B(col, col) + beta * C(col, col)*/
01593       Telement::ssmm(alpha, A(col,col), B(col, col),
01594                      beta_tmp, C(col,col));
01595       for (int ind = 0; ind < col; ind++)
01596         /* C(col, col) += alpha * A(ind, col)' * B(ind, col) */
01597         Telement::gemm(true, false, 
01598                        alpha, A(ind,col), B(ind, col),
01599                        1.0, C(col,col));
01600       for (int ind = col + 1; ind < A.ncols(); ind++)
01601         /* C(col, col) += alpha * A(col, ind) * B(col, ind)' */
01602         Telement::gemm(false, true,  
01603                        alpha, A(col, ind), B(col, ind),
01604                        1.0, C(col,col));
01605       /* Upper half */
01606       for (int row = 0; row < col; row++) {
01607         /* C(row, col) = alpha * A(row, row) * B(row, col) + 
01608          *               beta * C(row, col)
01609          */
01610         Telement::symm('L', 'U', 
01611                        alpha, A(row, row), B(row, col), 
01612                        beta_tmp, C(row, col));
01613         /* C(row, col) += alpha * A(row, col) * B(col, col)      */
01614         Telement::symm('R', 'U', 
01615                        alpha, B(col, col), A(row, col), 
01616                        1.0, C(row, col));
01617         for (int ind = 0; ind < row; ind++)
01618           /* C(row, col) += alpha * A(ind, row)' * B(ind, col)   */
01619           Telement::gemm(true, false,  
01620                          alpha, A(ind, row), B(ind, col),
01621                          1.0, C(row,col));
01622         for (int ind = row + 1; ind < col; ind++)
01623           /* C(row, col) += alpha * A(row, ind) * B(ind, col)    */
01624           Telement::gemm(false, false,  
01625                          alpha, A(row, ind), B(ind, col),
01626                          1.0, C(row,col));
01627         for (int ind = col + 1; ind < A.ncols(); ind++)
01628           /* C(row, col) += alpha * A(row, ind) * B(col, ind)'   */
01629           Telement::gemm(false, true,  
01630                          alpha, A(row, ind), B(col, ind),
01631                          1.0, C(row,col));
01632       }
01633       /* Lower half */
01634       Telement tmpSubMat;
01635       for (int row = col + 1; row < C.nrows(); row++) {
01636         Telement::transpose(C(row, col), tmpSubMat);
01637         /* tmpSubMat = alpha * B(col, col) * A(col, row) + 
01638          *             beta * tmpSubMat
01639          */
01640         Telement::symm('L', 'U', 
01641                        alpha, B(col, col), A(col, row), 
01642                        beta_tmp, tmpSubMat);
01643         /* tmpSubMat += alpha * B(col, row) * A(row, row)        */
01644         Telement::symm('R', 'U', 
01645                        alpha, A(row, row), B(col, row), 
01646                        1.0, tmpSubMat);
01647         for (int ind = 0; ind < col; ind++)
01648           /* tmpSubMat += alpha * B(ind, col)' * A(ind, row)     */
01649           Telement::gemm(true, false,  
01650                          alpha, B(ind, col), A(ind, row),
01651                          1.0, tmpSubMat);
01652         for (int ind = col + 1; ind < row; ind++)
01653           /* tmpSubMat += alpha * B(col, ind) * A(ind, row)      */
01654           Telement::gemm(false, false,  
01655                          alpha, B(col, ind), A(ind, row),
01656                          1.0, tmpSubMat);
01657         for (int ind = row + 1; ind < B.nrows(); ind++)
01658           /* tmpSubMat += alpha * B(col, ind) * A(row, ind)'     */
01659           Telement::gemm(false, true,  
01660                          alpha, B(col, ind), A(row, ind),
01661                          1.0, tmpSubMat);
01662         Telement::transpose(tmpSubMat, C(row, col));
01663       }
01664       MAT_OMP_END;
01665     }
01666     MAT_OMP_FINALIZE;
01667   } /* end ssmm */
01668   
01669 
01670     /* C = alpha * A * B + beta * C where A and B are symmetric
01671      * and only the upper triangle of C is computed.
01672      */
01673   template<class Treal, class Telement>
01674     void Matrix<Treal, Telement>:: 
01675     ssmm_upper_tr_only(const Treal alpha, 
01676                        const Matrix<Treal, Telement>& A, 
01677                        const Matrix<Treal, Telement>& B, 
01678                        const Treal beta, 
01679                        Matrix<Treal, Telement>& C) {
01680     assert(!A.is_empty());
01681     assert(!B.is_empty());
01682     if (C.is_empty()) {
01683       assert(beta == 0);
01684       C.resetRows(A.rows);
01685       C.resetCols(B.cols);
01686     }
01687     if (A.ncols() != B.nrows() || 
01688         A.nrows() != C.nrows() ||
01689         B.ncols() != C.ncols() ||
01690         A.nrows() != A.ncols() ||
01691         B.nrows() != B.ncols()) {
01692       throw Failure("Matrix<class Treal, class Telement>::"
01693                     "ssmm_upper_tr_only(Treal, "
01694                     "Matrix<Treal, Telement>, "
01695                     "Matrix<Treal, Telement>, Treal, "
01696                     "Matrix<Treal, Telement>): "
01697                     "Incorrect matrixdimensions for ssmm_upper_tr_only "
01698                     "multiplication");
01699     }
01700     if ((A.is_zero() || B.is_zero() || alpha == 0) && 
01701         (C.is_zero() || beta == 0)) {
01702       C = 0;
01703       return;
01704     }
01705     Treal beta_tmp = beta;
01706     if (C.is_zero()) {
01707       C.allocate();
01708       beta_tmp = 0;
01709     }
01710     if (A.is_zero() || B.is_zero() || alpha == 0) {
01711       C *= beta;
01712       return;
01713     }
01714     MAT_OMP_INIT;
01715 #ifdef _OPENMP
01716 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic) 
01717 #endif
01718     for (int col = 0; col < C.ncols(); col++) {
01719       MAT_OMP_START;
01720       /* Diagonal */
01721       /* C(col, col) = alpha * A(col, col) * B(col, col) + beta * C(col, col)*/
01722       Telement::ssmm_upper_tr_only(alpha, A(col,col), B(col, col),
01723                                    beta_tmp, C(col,col));
01724       for (int ind = 0; ind < col; ind++)
01725         /* C(col, col) += alpha * A(ind, col)' * B(ind, col) */
01726         Telement::gemm_upper_tr_only(true, false, 
01727                                      alpha, A(ind,col), B(ind, col),
01728                                      1.0, C(col,col));
01729       for (int ind = col + 1; ind < A.ncols(); ind++)
01730         /* C(col, col) += alpha * A(col, ind) * B(col, ind)' */
01731         Telement::gemm_upper_tr_only(false, true,  
01732                                      alpha, A(col, ind), B(col, ind),
01733                                      1.0, C(col,col));
01734       /* Upper half */
01735       for (int row = 0; row < col; row++) {
01736         /* C(row, col) = alpha * A(row, row) * B(row, col) + 
01737          *               beta * C(row, col)
01738          */
01739         Telement::symm('L', 'U', 
01740                        alpha, A(row, row), B(row, col), 
01741                        beta_tmp, C(row, col));
01742         /* C(row, col) += alpha * A(row, col) * B(col, col)      */
01743         Telement::symm('R', 'U', 
01744                        alpha, B(col, col), A(row, col), 
01745                        1.0, C(row, col));
01746         for (int ind = 0; ind < row; ind++)
01747           /* C(row, col) += alpha * A(ind, row)' * B(ind, col)   */
01748           Telement::gemm(true, false,  
01749                          alpha, A(ind, row), B(ind, col),
01750                          1.0, C(row,col));
01751         for (int ind = row + 1; ind < col; ind++)
01752           /* C(row, col) += alpha * A(row, ind) * B(ind, col)    */
01753           Telement::gemm(false, false,  
01754                          alpha, A(row, ind), B(ind, col),
01755                          1.0, C(row,col));
01756         for (int ind = col + 1; ind < A.ncols(); ind++)
01757           /* C(row, col) += alpha * A(row, ind) * B(col, ind)'   */
01758           Telement::gemm(false, true,  
01759                          alpha, A(row, ind), B(col, ind),
01760                          1.0, C(row,col));
01761       }
01762       MAT_OMP_END;
01763     }
01764     MAT_OMP_FINALIZE;
01765   } /* end ssmm_upper_tr_only */
01766   
01767 
01768 
01769   template<class Treal, class Telement>
01770     void Matrix<Treal, Telement>:: 
01771     trmm(const char side, const char uplo, const bool tA, const Treal alpha, 
01772          const Matrix<Treal, Telement>& A, 
01773          Matrix<Treal, Telement>& B) {
01774     assert(!A.is_empty());
01775     assert(!B.is_empty());
01776     if (alpha != 0 && !A.is_zero() && !B.is_zero())
01777       if (((side == 'R' && B.ncols() == A.nrows()) || 
01778            (side == 'L' && A.ncols() == B.nrows())) &&
01779           A.nrows() == A.ncols())
01780         if (uplo == 'U')
01781           if (!tA) { 
01782             if (side == 'R') {
01783               /* Last column must be calculated first                        */
01784               for (int col = B.ncols() - 1; col >= 0; col--)
01785                 for (int row = 0; row < B.nrows(); row++) {
01786                   /* Use first the diagonal element in A                     */
01787                   /* Otherwise the faster call to trmm cannot be utilized    */
01788                   Telement::trmm(side, uplo, tA, alpha, 
01789                                  A(col, col), B(row,col));
01790                   /* And the rest:                                           */
01791                   for (int ind = 0; ind < col; ind++)
01792                     Telement::gemm(false, tA, alpha, 
01793                                    B(row,ind), A(ind, col), 
01794                                    1.0, B(row,col));
01795                 }
01796             } /* end if (side == 'R')*/
01797             else {
01798               assert(side == 'L');
01799               /* First row must be calculated first                          */
01800               for (int row = 0; row < B.nrows(); row++ )
01801                 for (int col = 0; col < B.ncols(); col++) {
01802                   Telement::trmm(side, uplo, tA, alpha, 
01803                                  A(row,row), B(row,col));
01804                   for (int ind = row + 1 ; ind < B.nrows(); ind++)
01805                     Telement::gemm(tA, false, alpha,
01806                                    A(row,ind), B(ind,col),
01807                                    1.0, B(row,col));
01808                 }
01809             } /* end else (side == 'L')*/
01810           } /* end if (tA == false) */
01811           else {
01812             assert(tA == true);
01813             if (side == 'R') {
01814               /* First column must be calculated first                       */
01815               for (int col = 0; col < B.ncols(); col++)
01816                 for (int row = 0; row < B.nrows(); row++) {
01817                   Telement::trmm(side, uplo, tA, alpha,
01818                                  A(col,col), B(row,col));
01819                   for (int ind = col + 1; ind < A.ncols(); ind++)
01820                     Telement::gemm(false, tA, alpha,
01821                                    B(row,ind), A(col,ind),
01822                                    1.0, B(row,col));
01823                 }
01824             } /* end if (side == 'R')*/
01825             else {
01826               assert(side == 'L');
01827               /* Last row must be calculated first                           */
01828               for (int row = B.nrows() - 1; row >= 0; row--)
01829                 for (int col = 0; col < B.ncols(); col++) {
01830                   Telement::trmm(side, uplo, tA, alpha,
01831                                  A(row,row), B(row,col));
01832                   for (int ind = 0; ind < row; ind++)
01833                     Telement::gemm(tA, false, alpha,
01834                                    A(ind,row), B(ind,col),
01835                                    1.0, B(row,col));
01836                 }
01837             } /* end else (side == 'L')*/
01838           } /* end else (tA == true)*/
01839         else
01840           throw Failure("Matrix<class Treal, class Telement>::"
01841                         "trmm not implemented for lower triangular matrices");
01842       else
01843         throw Failure("Matrix<class Treal, class Telement>::trmm"
01844                       ": Incorrect matrix dimensions for multiplication");
01845     else
01846       B = 0;
01847   }
01848 
01849 
01850   template<class Treal, class Telement>
01851     Treal Matrix<Treal, Telement>::frobSquared() const {
01852     assert(!this->is_empty());
01853     if (this->is_zero()) 
01854       return 0;
01855     else {
01856       Treal sum(0);
01857       for (int i = 0; i < this->nElements(); i++)
01858         sum += this->elements[i].frobSquared();
01859       return sum;
01860     }
01861   } 
01862   
01863   template<class Treal, class Telement>
01864     Treal Matrix<Treal, Telement>::
01865     syFrobSquared() const {
01866     assert(!this->is_empty());
01867     if (this->nrows() != this->ncols())
01868       throw Failure("Matrix<Treal, Telement>::syFrobSquared(): "
01869                     "Matrix must be have equal number of rows "
01870                     "and cols to be symmetric");
01871     Treal sum(0);
01872     if (!this->is_zero()) {
01873       for (int col = 1; col < this->ncols(); col++)
01874         for (int row = 0; row < col; row++)
01875           sum += 2 * (*this)(row, col).frobSquared();
01876       for (int rc = 0; rc < this->ncols(); rc++)
01877         sum += (*this)(rc, rc).syFrobSquared();
01878     }
01879     return sum;
01880   }
01881   
01882   template<class Treal, class Telement>
01883     Treal Matrix<Treal, Telement>::
01884     frobSquaredDiff
01885     (const Matrix<Treal, Telement>& A,
01886      const Matrix<Treal, Telement>& B) {
01887     assert(!A.is_empty());
01888     assert(!B.is_empty());
01889     if (A.nrows() != B.nrows() || A.ncols() != B.ncols()) 
01890       throw Failure("Matrix<Treal, Telement>::"
01891                     "frobSquaredDiff: Incorrect matrix dimensions");
01892     Treal sum(0);
01893     if (!A.is_zero() && !B.is_zero()) 
01894       for (int i = 0; i < A.nElements(); i++)
01895         sum += Telement::frobSquaredDiff(A.elements[i],B.elements[i]);
01896     else if (A.is_zero() && !B.is_zero()) 
01897       sum = B.frobSquared();
01898     else if (!A.is_zero() && B.is_zero())
01899       sum = A.frobSquared();
01900     /* sum is already zero if A.is_zero() && B.is_zero() */
01901     return sum;
01902   }  
01903   
01904   template<class Treal, class Telement>
01905     Treal Matrix<Treal, Telement>::
01906     syFrobSquaredDiff
01907     (const Matrix<Treal, Telement>& A,
01908      const Matrix<Treal, Telement>& B) {
01909     assert(!A.is_empty());
01910     assert(!B.is_empty());
01911     if (A.nrows() != B.nrows() ||
01912         A.ncols() != B.ncols() || 
01913         A.nrows() != A.ncols()) 
01914       throw Failure("Matrix<Treal, Telement>::syFrobSquaredDiff: "
01915                     "Incorrect matrix dimensions");
01916     Treal sum(0);
01917     if (!A.is_zero() && !B.is_zero()) {
01918       for (int col = 1; col < A.ncols(); col++)
01919         for (int row = 0; row < col; row++)
01920           sum += 2 * Telement::frobSquaredDiff(A(row, col), B(row, col));
01921       for (int rc = 0; rc < A.ncols(); rc++)
01922         sum += Telement::syFrobSquaredDiff(A(rc, rc),B(rc, rc));
01923     }
01924     else if (A.is_zero() && !B.is_zero()) 
01925       sum = B.syFrobSquared();
01926     else if (!A.is_zero() && B.is_zero())
01927       sum = A.syFrobSquared();
01928     /* sum is already zero if A.is_zero() && B.is_zero() */
01929     return sum;
01930   }
01931   
01932   template<class Treal, class Telement> 
01933     Treal Matrix<Treal, Telement>::
01934     trace() const {
01935     assert(!this->is_empty());
01936     if (this->nrows() != this->ncols())
01937       throw Failure("Matrix<Treal, Telement>::trace(): "
01938                     "Matrix must be quadratic");          
01939     if (this->is_zero())
01940       return 0;
01941     else {
01942       Treal sum = 0;
01943       for (int rc = 0; rc < this->ncols(); rc++)
01944         sum += (*this)(rc,rc).trace();
01945       return sum;
01946     }
01947   }
01948   
01949   template<class Treal, class Telement> 
01950     Treal Matrix<Treal, Telement>::
01951     trace_ab(const Matrix<Treal, Telement>& A,
01952              const Matrix<Treal, Telement>& B) {
01953     assert(!A.is_empty());
01954     assert(!B.is_empty());
01955     if (A.nrows() != B.ncols() || A.ncols() != B.nrows()) 
01956       throw Failure("Matrix<Treal, Telement>::trace_ab: "
01957                     "Wrong matrix dimensions for trace(A * B)"); 
01958     Treal tr = 0;
01959     if (!A.is_zero() && !B.is_zero())
01960       for (int rc = 0; rc < A.nrows(); rc++)
01961         for (int colA = 0; colA < A.ncols(); colA++)
01962           tr += Telement::trace_ab(A(rc,colA), B(colA, rc));
01963     return tr;
01964   } 
01965   
01966   template<class Treal, class Telement> 
01967     Treal  Matrix<Treal, Telement>::
01968     trace_aTb(const Matrix<Treal, Telement>& A,
01969               const Matrix<Treal, Telement>& B) {
01970     assert(!A.is_empty());
01971     assert(!B.is_empty());
01972     if (A.ncols() != B.ncols() || A.nrows() != B.nrows())
01973       throw Failure("Matrix<Treal, Telement>::trace_aTb: "
01974                     "Wrong matrix dimensions for trace(A' * B)"); 
01975     Treal tr = 0;
01976     if (!A.is_zero() && !B.is_zero())
01977       for (int rc = 0; rc < A.ncols(); rc++)
01978         for (int rowB = 0; rowB < B.nrows(); rowB++)
01979           tr += Telement::trace_aTb(A(rowB,rc), B(rowB, rc));
01980     return tr;
01981   }
01982   
01983   template<class Treal, class Telement> 
01984     Treal Matrix<Treal, Telement>::
01985     sy_trace_ab(const Matrix<Treal, Telement>& A,
01986                 const Matrix<Treal, Telement>& B) {
01987     assert(!A.is_empty());
01988     assert(!B.is_empty());
01989     if (A.nrows() != B.ncols() || A.ncols() != B.nrows() || 
01990         A.nrows() != A.ncols()) 
01991       throw Failure("Matrix<Treal, Telement>::sy_trace_ab: "
01992                     "Wrong matrix dimensions for symmetric trace(A * B)");
01993     Treal tr = 0;
01994     if (!A.is_zero() && !B.is_zero()) {
01995       /* Diagonal first */
01996       for (int rc = 0; rc < A.nrows(); rc++)
01997         tr += Telement::sy_trace_ab(A(rc,rc), B(rc, rc));
01998       /* Using that trace of transpose is equal to that without transpose: */
01999       for (int colA = 1; colA < A.ncols(); colA++)
02000         for (int rowA = 0; rowA < colA; rowA++)
02001           tr += 2 * Telement::trace_aTb(A(rowA, colA), B(rowA, colA));
02002     }
02003     return tr;
02004   }
02005 
02006   template<class Treal, class Telement>
02007     void Matrix<Treal, Telement>:: 
02008     add(const Treal alpha,   /* B += alpha * A */
02009         const Matrix<Treal, Telement>& A, 
02010         Matrix<Treal, Telement>& B) {
02011     assert(!A.is_empty());
02012     assert(!B.is_empty());
02013     if (A.nrows() != B.nrows() || A.ncols() != B.ncols()) 
02014       throw Failure("Matrix<Treal, Telement>::add: "
02015                     "Wrong matrix dimensions for addition");
02016     if (!A.is_zero() && alpha != 0) {
02017       if (B.is_zero())
02018         B.allocate();
02019       for (int ind = 0; ind < A.nElements(); ind++)
02020         Telement::add(alpha, A.elements[ind], B.elements[ind]);
02021     }
02022   }
02023 
02024   template<class Treal, class Telement>
02025     void Matrix<Treal, Telement>:: 
02026     assign(Treal const  alpha,   /* *this = alpha * A */
02027            Matrix<Treal, Telement> const & A) {
02028     assert(!A.is_empty());
02029     if (this->is_empty()) {
02030       this->resetRows(A.rows);
02031       this->resetCols(A.cols);
02032     }
02033     *this = 0;
02034     Matrix<Treal, Telement>:: 
02035       add(alpha, A, *this);
02036   }
02037 
02038 
02039   /********** Help functions for thresholding */
02040 
02041   
02042   template<class Treal, class Telement>
02043     void Matrix<Treal, Telement>::
02044     getFrobSqLowestLevel(std::vector<Treal> & frobsq) const {
02045     if (!this->is_zero())
02046       for (int ind = 0; ind < this->nElements(); ind++)
02047         this->elements[ind].getFrobSqLowestLevel(frobsq);
02048   }
02049 
02050   template<class Treal, class Telement>
02051     void Matrix<Treal, Telement>::
02052     frobThreshLowestLevel
02053     (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix) {
02054     assert(!this->is_empty());
02055     bool thisMatIsZero = true;
02056     if (ErrorMatrix) {
02057       assert(!ErrorMatrix->is_empty());
02058       bool EMatIsZero = true;
02059       if (!ErrorMatrix->is_zero() || !this->is_zero()) {
02060         if (ErrorMatrix->is_zero())
02061           ErrorMatrix->allocate();
02062         if (this->is_zero())
02063           this->allocate();
02064         for (int ind = 0; ind < this->nElements(); ind++) {
02065           this->elements[ind].
02066             frobThreshLowestLevel(threshold, &ErrorMatrix->elements[ind]);
02067           thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
02068           EMatIsZero    = EMatIsZero && ErrorMatrix->elements[ind].is_zero();
02069         }
02070         if (thisMatIsZero)
02071           this->clear();
02072         if (EMatIsZero)
02073           ErrorMatrix->clear();
02074       }
02075     }      
02076     else
02077       if (!this->is_zero()) {
02078         for (int ind = 0; ind < this->nElements(); ind++) {
02079           this->elements[ind].frobThreshLowestLevel(threshold, 0);
02080           thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
02081         }
02082         if (thisMatIsZero)
02083           this->clear();
02084       }
02085   }
02086 
02087   template<class Treal, class Telement>
02088     void Matrix<Treal, Telement>::
02089     getFrobSqElementLevel(std::vector<Treal> & frobsq) const {
02090     if (!this->is_zero())
02091       for (int ind = 0; ind < this->nElements(); ind++)
02092         this->elements[ind].getFrobSqElementLevel(frobsq);
02093   }
02094 
02095   template<class Treal, class Telement>
02096     void Matrix<Treal, Telement>::
02097     frobThreshElementLevel
02098     (Treal const threshold, Matrix<Treal, Telement> * ErrorMatrix) {
02099     assert(!this->is_empty());
02100     bool thisMatIsZero = true;
02101     if (ErrorMatrix) {
02102       assert(!ErrorMatrix->is_empty());
02103       bool EMatIsZero = true;
02104       if (!ErrorMatrix->is_zero() || !this->is_zero()) {
02105         if (ErrorMatrix->is_zero())
02106           ErrorMatrix->allocate();
02107         if (this->is_zero())
02108           this->allocate();
02109         for (int ind = 0; ind < this->nElements(); ind++) {
02110           this->elements[ind].
02111             frobThreshElementLevel(threshold, &ErrorMatrix->elements[ind]);
02112           thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
02113           EMatIsZero    = EMatIsZero && ErrorMatrix->elements[ind].is_zero();
02114         }
02115         if (thisMatIsZero)
02116           this->clear();
02117         if (EMatIsZero)
02118           ErrorMatrix->clear();
02119       }
02120     }      
02121     else
02122       if (!this->is_zero()) {
02123         for (int ind = 0; ind < this->nElements(); ind++) {
02124           this->elements[ind].frobThreshElementLevel(threshold, 0);
02125           thisMatIsZero = thisMatIsZero && this->elements[ind].is_zero();
02126         }
02127         if (thisMatIsZero)
02128           this->clear();
02129       }    
02130   }
02131   
02132 
02133 
02134   template<class Treal, class Telement>
02135     void Matrix<Treal, Telement>::assignFrobNormsLowestLevel
02136     ( Matrix<Treal, Matrix<Treal, Telement> > const & A ) {
02137     if (!A.is_zero()) {
02138       if ( this->is_zero() )
02139         this->allocate();
02140       assert( this->nElements() == A.nElements() );
02141       for (int ind = 0; ind < this->nElements(); ind++)
02142         this->elements[ind].assignFrobNormsLowestLevel(A[ind]);
02143     }
02144     else
02145       this->clear();
02146   } 
02147 
02148   template<class Treal, class Telement>
02149     void Matrix<Treal, Telement>::syAssignFrobNormsLowestLevel
02150     ( Matrix<Treal, Matrix<Treal, Telement> > const & A ) {
02151     assert(!A.is_empty());
02152     if (A.nrows() != A.ncols())
02153       throw Failure("Matrix<Treal, Telement>::syAssignFrobNormsLowestLevel(...): "
02154                     "Matrix must be have equal number of rows "
02155                     "and cols to be symmetric");
02156     if (!A.is_zero()) {
02157       if ( this->is_zero() )
02158         this->allocate();
02159       assert( this->nElements() == A.nElements() );
02160       for (int col = 1; col < this->ncols(); col++)
02161         for (int row = 0; row < col; row++)
02162           (*this)(row, col).assignFrobNormsLowestLevel( A(row,col) );
02163       for (int rc = 0; rc < this->ncols(); rc++)
02164         (*this)(rc, rc).syAssignFrobNormsLowestLevel( A(rc,rc) );
02165     }
02166     else 
02167       this->clear();    
02168   }
02169 
02170   template<class Treal, class Telement>
02171     void Matrix<Treal, Telement>::assignDiffFrobNormsLowestLevel
02172     ( Matrix<Treal, Matrix<Treal, Telement> > const & A,
02173       Matrix<Treal, Matrix<Treal, Telement> > const & B ) {
02174     if ( A.is_zero() && B.is_zero() ) {
02175       // Both A and B are zero
02176       this->clear();
02177       return;
02178     }
02179     // At least one of A and B is nonzero
02180     if ( this->is_zero() )
02181       this->allocate();
02182     if ( !A.is_zero() && !B.is_zero() ) {
02183       // Both are nonzero
02184       assert( this->nElements() == A.nElements() );
02185       assert( this->nElements() == B.nElements() );
02186       for (int ind = 0; ind < this->nElements(); ind++)
02187         this->elements[ind].assignDiffFrobNormsLowestLevel( A[ind], B[ind] );
02188       return;
02189     }
02190     if ( !A.is_zero() ) {
02191       // A is nonzero
02192       this->assignFrobNormsLowestLevel( A );
02193       return;
02194     }
02195     if ( !B.is_zero() ) {
02196       // B is nonzero
02197       this->assignFrobNormsLowestLevel( B );
02198       return;
02199     }
02200   }
02201   template<class Treal, class Telement>
02202     void Matrix<Treal, Telement>::syAssignDiffFrobNormsLowestLevel
02203     ( Matrix<Treal, Matrix<Treal, Telement> > const & A,
02204       Matrix<Treal, Matrix<Treal, Telement> > const & B ) {
02205     if ( A.is_zero() && B.is_zero() ) {
02206       // Both A and B are zero
02207       this->clear();
02208       return;
02209     }
02210     // At least one of A and B is nonzero
02211     if ( this->is_zero() )
02212       this->allocate();
02213     if ( !A.is_zero() && !B.is_zero() ) {
02214       // Both are nonzero
02215       assert( this->nElements() == A.nElements() );
02216       assert( this->nElements() == B.nElements() );
02217       for (int col = 1; col < this->ncols(); col++)
02218         for (int row = 0; row < col; row++)
02219           (*this)(row, col).assignDiffFrobNormsLowestLevel( A(row,col), B(row,col) );
02220       for (int rc = 0; rc < this->ncols(); rc++)
02221         (*this)(rc, rc).syAssignDiffFrobNormsLowestLevel( A(rc,rc), B(rc,rc) );
02222       return;
02223     }
02224     if ( !A.is_zero() ) {
02225       // A is nonzero
02226       this->syAssignFrobNormsLowestLevel( A );
02227       return;
02228     }
02229     if ( !B.is_zero() ) {
02230       // B is nonzero
02231       this->syAssignFrobNormsLowestLevel( B );
02232       return;
02233     }
02234   } 
02235 
02236 
02237 
02238   template<class Treal, class Telement>
02239     void Matrix<Treal, Telement>::
02240     truncateAccordingToSparsityPattern( Matrix<Treal, Matrix<Treal, Telement> > & A ) const {
02241     if ( this->is_zero() ) {
02242       A.clear();
02243     }
02244     else {
02245       assert( !A.is_zero() );
02246       assert( this->nElements() == A.nElements() );
02247       for (int ind = 0; ind < this->nElements(); ind++)
02248         this->elements[ind].truncateAccordingToSparsityPattern( A[ind] );
02249     }
02250   }  
02251   
02252 
02253 
02254   /********** End of help functions for thresholding */
02255   
02256   /* C = beta * C + alpha * A * B where only the upper triangle of C is */
02257   /* referenced and updated */
02258   template<class Treal, class Telement>
02259     void Matrix<Treal, Telement>::
02260     gemm_upper_tr_only(const bool tA, const bool tB, const Treal alpha, 
02261                        const Matrix<Treal, Telement>& A, 
02262                        const Matrix<Treal, Telement>& B, 
02263                        const Treal beta, 
02264                        Matrix<Treal, Telement>& C) {
02265     assert(!A.is_empty());
02266     assert(!B.is_empty());
02267     if (C.is_empty()) {
02268       assert(beta == 0);
02269       if (tA)
02270         C.resetRows(A.cols);
02271       else
02272         C.resetRows(A.rows);
02273       if (tB)
02274         C.resetCols(B.rows);
02275       else
02276         C.resetCols(B.cols);
02277     }      
02278     if ((A.is_zero() || B.is_zero() || alpha == 0) && 
02279         (C.is_zero() || beta == 0))
02280       C = 0;
02281     else {
02282       Treal beta_tmp = beta;
02283       if (C.is_zero()) {
02284         C.allocate();
02285         beta_tmp = 0;
02286       }
02287       if (!A.is_zero() && !B.is_zero() && alpha != 0) {
02288         if (!tA && !tB) 
02289           if (A.ncols() == B.nrows() && 
02290               A.nrows() == C.nrows() &&
02291               B.ncols() == C.ncols()) {
02292             for (int col = 0; col < C.ncols(); col++) {
02293               Telement::gemm_upper_tr_only(tA, tB, alpha, 
02294                                            A(col, 0), B(0, col), 
02295                                            beta_tmp, 
02296                                            C(col, col));
02297               for(int cola = 1; cola < A.ncols(); cola++)
02298                 Telement::gemm_upper_tr_only(tA, tB, alpha, 
02299                                              A(col, cola), B(cola, col), 
02300                                              1.0, 
02301                                              C(col,col));
02302               for (int row = 0; row < col; row++) {
02303                 Telement::gemm(tA, tB, alpha, 
02304                                A(row, 0), B(0, col), 
02305                                beta_tmp, 
02306                                C(row,col));
02307                 for(int cola = 1; cola < A.ncols(); cola++)
02308                   Telement::gemm(tA, tB, alpha, 
02309                                  A(row, cola), B(cola, col), 
02310                                  1.0, 
02311                                  C(row,col));
02312               }
02313             }
02314           }
02315           else 
02316             throw Failure("Matrix<class Treal, class Telement>::"
02317                           "gemm_upper_tr_only(bool, bool, Treal, "
02318                           "Matrix<Treal, Telement>, "
02319                           "Matrix<Treal, Telement>, "
02320                           "Treal, Matrix<Treal, Telement>): "
02321                           "Incorrect matrixdimensions for multiplication");
02322         else if (tA && !tB)
02323           if (A.nrows() == B.nrows() &&
02324               A.ncols() == C.nrows() &&
02325               B.ncols() == C.ncols()) {
02326             for (int col = 0; col < C.ncols(); col++) {
02327               Telement::gemm_upper_tr_only(tA, tB, alpha,
02328                                            A(0,col), B(0,col),
02329                                            beta_tmp,
02330                                            C(col,col));
02331               for(int cola = 1; cola < A.nrows(); cola++)
02332                 Telement::gemm_upper_tr_only(tA, tB, alpha, 
02333                                              A(cola, col), B(cola, col), 
02334                                              1.0, 
02335                                              C(col, col));
02336               for (int row = 0; row < col; row++) {
02337                 Telement::gemm(tA, tB, alpha,
02338                                A(0,row), B(0,col),
02339                                beta_tmp,
02340                                C(row,col));
02341                 for(int cola = 1; cola < A.nrows(); cola++)
02342                   Telement::gemm(tA, tB, alpha, 
02343                                  A(cola, row), B(cola, col), 
02344                                  1.0, 
02345                                  C(row,col));
02346               }
02347             }
02348           }
02349           else
02350             throw Failure("Matrix<class Treal, class Telement>::"
02351                           "gemm_upper_tr_only(bool, bool, Treal, "
02352                           "Matrix<Treal, Telement>, "
02353                           "Matrix<Treal, Telement>, "
02354                           "Treal, Matrix<Treal, Telement>): "
02355                           "Incorrect matrixdimensions for multiplication");
02356         else if (!tA && tB)
02357           if (A.ncols() == B.ncols() && 
02358               A.nrows() == C.nrows() &&
02359               B.nrows() == C.ncols()) {
02360             for (int col = 0; col < C.ncols(); col++) {
02361               Telement::gemm_upper_tr_only(tA, tB, alpha, 
02362                                            A(col, 0), B(col, 0), 
02363                                            beta_tmp, 
02364                                            C(col,col));
02365               for(int cola = 1; cola < A.ncols(); cola++)
02366                 Telement::gemm_upper_tr_only(tA, tB, alpha, 
02367                                              A(col, cola), B(col, cola), 
02368                                              1.0, 
02369                                              C(col,col));
02370               for (int row = 0; row < col; row++) {
02371                 Telement::gemm(tA, tB, alpha, 
02372                                A(row, 0), B(col, 0), 
02373                                beta_tmp, 
02374                                C(row,col));
02375                 for(int cola = 1; cola < A.ncols(); cola++)
02376                   Telement::gemm(tA, tB, alpha, 
02377                                  A(row, cola), B(col, cola), 
02378                                  1.0, 
02379                                  C(row,col));
02380               }
02381             }
02382           }
02383           else 
02384             throw Failure("Matrix<class Treal, class Telement>::"
02385                           "gemm_upper_tr_only(bool, bool, Treal, "
02386                           "Matrix<Treal, Telement>, "
02387                           "Matrix<Treal, Telement>, "
02388                           "Treal, Matrix<Treal, Telement>): "
02389                           "Incorrect matrixdimensions for multiplication");
02390         else if (tA && tB)
02391           if (A.nrows() == B.ncols() && 
02392               A.ncols() == C.nrows() &&
02393               B.nrows() == C.ncols()) {
02394             for (int col = 0; col < C.ncols(); col++) {
02395               Telement::gemm_upper_tr_only(tA, tB, alpha, 
02396                                            A(0, col), B(col, 0), 
02397                                            beta_tmp, 
02398                                            C(col,col));
02399               for(int cola = 1; cola < A.nrows(); cola++)
02400                 Telement::gemm_upper_tr_only(tA, tB, alpha, 
02401                                              A(cola, col), B(col, cola), 
02402                                              1.0, 
02403                                              C(col,col));
02404               for (int row = 0; row < col; row++) {
02405                 Telement::gemm(tA, tB, alpha, 
02406                                A(0, row), B(col, 0), 
02407                                beta_tmp, 
02408                                C(row,col));
02409                 for(int cola = 1; cola < A.nrows(); cola++)
02410                   Telement::gemm(tA, tB, alpha, 
02411                                  A(cola, row), B(col, cola), 
02412                                  1.0, 
02413                                  C(row,col));
02414               }
02415             }
02416           }
02417           else 
02418             throw Failure("Matrix<class Treal, class Telement>::"
02419                           "gemm_upper_tr_only(bool, bool, Treal, "
02420                           "Matrix<Treal, Telement>, "
02421                           "Matrix<Treal, Telement>, Treal, "
02422                           "Matrix<Treal, Telement>): "
02423                           "Incorrect matrixdimensions for multiplication");
02424         else throw Failure("Matrix<class Treal, class Telement>::"
02425                            "gemm_upper_tr_only(bool, bool, Treal, "
02426                            "Matrix<Treal, Telement>, "
02427                            "Matrix<Treal, Telement>, Treal, "
02428                            "Matrix<Treal, Telement>):"
02429                            "Very strange error!!");
02430       }    
02431       else
02432         C *= beta;
02433     }
02434   }
02435 
02436   /* A = alpha * A * Z or A = alpha * Z * A where A is symmetric, */
02437   /* Z is upper triangular and */
02438   /* only the upper triangle of the result is calculated */
02439   /* side == 'R' for A = alpha * A * Z */
02440   /* side == 'L' for A = alpha * Z * A */
02441   template<class Treal, class Telement>
02442     void Matrix<Treal, Telement>:: 
02443     sytr_upper_tr_only(char const side, const Treal alpha,
02444                        Matrix<Treal, Telement>& A,
02445                        const Matrix<Treal, Telement>& Z) {
02446     assert(!A.is_empty());
02447     assert(!Z.is_empty());
02448     if (alpha != 0 && !A.is_zero() && !Z.is_zero()) {
02449       if (A.nrows() == A.ncols() && 
02450           Z.nrows() == Z.ncols() && 
02451           A.nrows() == Z.nrows()) {
02452         if (side == 'R') {
02453           /* Last column must be calculated first                       */
02454           for (int col = A.ncols() - 1; col >= 0; col--) {
02455             // A(col, col) = alpha * A(col, col) * Z(col, col)
02456             Telement::sytr_upper_tr_only(side, alpha, 
02457                                          A(col, col), Z(col, col));
02458             for (int ind = 0; ind < col; ind++) {
02459               // A(col, col) += alpha * A(ind, col)' * Z(ind, col)
02460               Telement::gemm_upper_tr_only(true, false, alpha, A(ind, col), 
02461                                            Z(ind, col), 1.0, A(col, col));
02462             }
02463             /* Last row must be calculated first                       */
02464             for (int row = col - 1; row >= 0; row--) {
02465               // A(row, col) = alpha * A(row, col) * Z(col, col);
02466               Telement::trmm(side, 'U', false, 
02467                              alpha, Z(col, col), A(row, col));
02468               // A(row, col) += alpha * A(row, row) * Z(row, col);
02469               Telement::symm('L', 'U', alpha, A(row, row), Z(row, col), 
02470                              1.0, A(row, col));
02471               for (int ind = 0; ind < row; ind++) {
02472                 // A(row, col) += alpha * A(ind, row)' * Z(ind, col);
02473                 Telement::gemm(true, false, alpha, A(ind, row), Z(ind, col),
02474                                1.0, A(row, col));
02475               }
02476               for (int ind = row + 1; ind < col; ind++) {
02477                 // A(row, col) += alpha * A(row, ind) * Z(ind, col);
02478                 Telement::gemm(false, false, alpha, A(row, ind), Z(ind, col),
02479                                1.0, A(row, col));
02480               }
02481             }
02482           }
02483         }
02484         else { /* side == 'L' */
02485           assert(side == 'L');
02486           /* First column must be calculated first                       */
02487           for (int col = 0; col < A.ncols(); col++) {
02488             /* First row must be calculated first                       */
02489             for (int row = 0; row < col; row++) {
02490               // A(row, col) = alpha * Z(row, row) * A(row, col)
02491               Telement::trmm(side, 'U', false, alpha, 
02492                              Z(row, row), A(row, col));
02493               // A(row, col) += alpha * Z(row, col) * A(col, col)
02494               Telement::symm('R', 'U', alpha, A(col, col), Z(row, col), 
02495                              1.0, A(row, col));
02496               for (int ind = row + 1; ind < col; ind++)
02497                 // A(row, col) += alpha * Z(row, ind) * A(ind, col)
02498                 Telement::gemm(false, false, alpha, Z(row, ind), A(ind, col),
02499                                1.0, A(row, col));
02500               for (int ind = col + 1; ind < A.nrows(); ind++)
02501                 // A(row, col) += alpha * Z(row, ind) * A(col, ind)'
02502                 Telement::gemm(false, true, alpha, Z(row, ind), A(col, ind),
02503                                1.0, A(row, col));
02504             }
02505             Telement::sytr_upper_tr_only(side, alpha, 
02506                                          A(col, col), Z(col, col));
02507             for (int ind = col + 1; ind < A.ncols(); ind++)
02508               Telement::gemm_upper_tr_only(false, true, alpha, Z(col, ind), 
02509                                            A(col, ind), 1.0, A(col, col));
02510           }
02511         }
02512       }
02513       else      
02514         throw Failure("Matrix<class Treal, class Telement>::"
02515                       "sytr_upper_tr_only: Incorrect matrix dimensions "
02516                       "for symmetric-triangular multiplication");
02517     }     
02518     else
02519       A = 0;
02520   }
02521 
02522   /* The result B is assumed to be symmetric, i.e. only upper triangle */
02523   /* calculated and hence only upper triangle of input B is used. */
02524   template<class Treal, class Telement>
02525     void Matrix<Treal, Telement>:: 
02526     trmm_upper_tr_only(const char side, const char uplo, 
02527                        const bool tA, const Treal alpha, 
02528                        const Matrix<Treal, Telement>& A, 
02529                        Matrix<Treal, Telement>& B) {
02530     assert(!A.is_empty());
02531     assert(!B.is_empty());
02532     if (alpha != 0 && !A.is_zero() && !B.is_zero())
02533       if (((side == 'R' && B.ncols() == A.nrows()) || 
02534            (side == 'L' && A.ncols() == B.nrows())) &&
02535           A.nrows() == A.ncols())
02536         if (uplo == 'U')
02537           if (!tA) { 
02538             throw Failure("Matrix<Treal, class Telement>::"
02539                           "trmm_upper_tr_only : "
02540                           "not possible for upper triangular not transposed "
02541                           "matrices A since lower triangle of B is needed");
02542           } /* end if (tA == false) */
02543           else {
02544             assert(tA == true);
02545             if (side == 'R') {
02546               /* First column must be calculated first                       */
02547               for (int col = 0; col < B.ncols(); col++) {
02548                 Telement::trmm_upper_tr_only(side, uplo, tA, alpha,
02549                                              A(col,col), B(col,col));
02550                 for (int ind = col + 1; ind < A.ncols(); ind++)
02551                   Telement::gemm_upper_tr_only(false, tA, alpha,
02552                                                B(col,ind), A(col,ind),
02553                                                1.0, B(col,col));
02554                 for (int row = 0; row < col; row++) {
02555                   Telement::trmm(side, uplo, tA, alpha,
02556                                  A(col,col), B(row,col));
02557                   for (int ind = col + 1; ind < A.ncols(); ind++)
02558                     Telement::gemm(false, tA, alpha,
02559                                    B(row,ind), A(col,ind),
02560                                    1.0, B(row,col));
02561                 }
02562               }
02563             } /* end if (side == 'R')*/
02564             else {
02565               assert(side == 'L');
02566               /* Last row must be calculated first                           */
02567               for (int row = B.nrows() - 1; row >= 0; row--) {
02568                 Telement::trmm_upper_tr_only(side, uplo, tA, alpha,
02569                                              A(row,row), B(row,row));
02570                 for (int ind = 0; ind < row; ind++)
02571                   Telement::gemm_upper_tr_only(tA, false, alpha,
02572                                                A(ind,row), B(ind,row),
02573                                                1.0, B(row,row));
02574                 for (int col = row + 1; col < B.ncols(); col++) {
02575                   Telement::trmm(side, uplo, tA, alpha,
02576                                  A(row,row), B(row,col));
02577                   for (int ind = 0; ind < row; ind++)
02578                     Telement::gemm(tA, false, alpha,
02579                                    A(ind,row), B(ind,col),
02580                                    1.0, B(row,col));
02581                 }
02582               }
02583             } /* end else (side == 'L')*/
02584           } /* end else (tA == true)*/
02585         else
02586           throw Failure("Matrix<class Treal, class Telement>::"
02587                         "trmm_upper_tr_only not implemented for lower "
02588                         "triangular matrices");
02589       else
02590         throw Failure("Matrix<class Treal, class Telement>::"
02591                       "trmm_upper_tr_only: Incorrect matrix dimensions "
02592                       "for multiplication");
02593     else
02594       B = 0;
02595   }
02596 
02597   /* A = Z' * A * Z or A = Z * A * Z' */
02598   /* where A is symmetric and Z is (nonunit) upper triangular */
02599   /* side == 'R' for A = Z' * A * Z */
02600   /* side == 'L' for A = Z * A * Z' */
02601   template<class Treal, class Telement>
02602     void Matrix<Treal, Telement>:: 
02603     trsytriplemm(char const side, 
02604                  const Matrix<Treal, Telement>& Z,
02605                  Matrix<Treal, Telement>& A) {
02606     if (side == 'R') {
02607       Matrix<Treal, Telement>::
02608 	sytr_upper_tr_only('R', 1.0, A, Z);
02609       Matrix<Treal, Telement>::
02610 	trmm_upper_tr_only('L', 'U', true, 1.0, Z, A);
02611     }
02612     else {
02613       assert(side == 'L');
02614       Matrix<Treal, Telement>::
02615 	sytr_upper_tr_only('L', 1.0, A, Z);
02616       Matrix<Treal, Telement>::
02617 	trmm_upper_tr_only('R', 'U', true, 1.0, Z, A);
02618     }
02619   }
02620 
02621   template<class Treal, class Telement>
02622     Treal Matrix<Treal, Telement>::
02623     frob_squared_thresh(Treal const threshold,
02624                         Matrix<Treal, Telement> * ErrorMatrix) {
02625     assert(!this->is_empty());
02626     if (ErrorMatrix && ErrorMatrix->is_empty()) {
02627       ErrorMatrix->resetRows(this->rows);
02628       ErrorMatrix->resetCols(this->cols);
02629     }
02630     assert(threshold >= (Treal)0.0);
02631     if (threshold == (Treal)0.0)
02632       return 0;
02633     
02634     std::vector<Treal> frobsq_norms;
02635     getFrobSqLowestLevel(frobsq_norms);
02636     /* Sort in ascending order */
02637     std::sort(frobsq_norms.begin(), frobsq_norms.end());
02638     int lastRemoved = -1;
02639     Treal frobsqSum = 0;
02640     int nnzBlocks = frobsq_norms.size();
02641     while(lastRemoved < nnzBlocks - 1 && frobsqSum < threshold) {
02642       ++lastRemoved;
02643       frobsqSum += frobsq_norms[lastRemoved];
02644     }
02645 
02646     /* Check if entire matrix will be removed */
02647     if (lastRemoved == nnzBlocks - 1 && frobsqSum < threshold) {
02648       if (ErrorMatrix) 
02649         Matrix<Treal, Telement>::swap(*this, *ErrorMatrix);
02650       else
02651         *this = 0;
02652     }
02653     else {
02654       frobsqSum -= frobsq_norms[lastRemoved];
02655       --lastRemoved;
02656       while(lastRemoved > -1 && frobsq_norms[lastRemoved] == 
02657             frobsq_norms[lastRemoved + 1]) {
02658         frobsqSum -= frobsq_norms[lastRemoved];
02659         --lastRemoved;
02660       }
02661       if (lastRemoved > -1) {
02662         Treal threshLowestLevel = 
02663           (frobsq_norms[lastRemoved + 1] + 
02664            frobsq_norms[lastRemoved]) / 2;
02665         this->frobThreshLowestLevel(threshLowestLevel, ErrorMatrix);
02666       }
02667     }
02668     return frobsqSum;
02669   }
02670 
02671   template<class Treal, class Telement> 
02672     void Matrix<Treal, Telement>::
02673     syInch(const Matrix<Treal, Telement>& A,
02674            Matrix<Treal, Telement>& Z,
02675            const Treal threshold, const side looking,
02676            const inchversion version) {
02677     assert(!A.is_empty());
02678     if (A.nrows() != A.ncols()) 
02679       throw Failure("Matrix<Treal, Telement>::sy_inch: "
02680                     "Matrix must be quadratic!");
02681     if (A.is_zero())
02682       throw Failure("Matrix<Treal>::sy_inch: Matrix is zero! "
02683                     "Not possible to compute inverse cholesky.");
02684     if (version == stable) /* STABILIZED */
02685       throw Failure("Matrix<Treal>::sy_inch: Only unstable "
02686                     "version of sy_inch implemented.");
02687     Treal myThresh = 0;
02688     if (threshold != 0)
02689       myThresh = threshold / (Z.ncols() * Z.nrows());
02690     Z.resetRows(A.rows);
02691     Z.resetCols(A.cols);
02692     Z.allocate();
02693     
02694     if (looking == left) {  /* LEFT-LOOKING INCH */
02695       if (threshold != 0)
02696         throw Failure("Matrix<Treal, Telement>::syInch: "
02697                       "Thresholding not implemented for left-looking inch.");
02698       /* Left and unstable */
02699       Telement::syInch(A(0,0), Z(0,0), threshold, looking, version);
02700       Telement Ptmp;//, tmp;
02701       for (int i = 1; i < Z.ncols(); i++) {
02702         for (int j = 0; j < i; j++) {   
02703           /* Z(k,i) is nonzero for k = 0, 1, ...,j - 2, j - 1, i         */
02704           /* and Z(i,i) = I                       (yes it is i ^)        */
02705           Ptmp = A(j,i);                    /* (Z(i,i) == I)             */
02706           for (int k = 0; k < j; k++)       /* Ptmp = A(j,:) * Z(:,i)  */
02707             Telement::gemm(true, false, 1.0, /* SYMMETRY USED */
02708                            A(k,j), Z(k,i), 1.0, Ptmp);
02709           Telement::trmm('L', 'U', true, 1.0, Z(j,j), Ptmp);
02710         
02711           for (int k = 0; k < j; k++)       /* Z(:,i) -= Z(:,j) * Ptmp   */
02712             Telement::gemm(false, false, -1.0,
02713                            Z(k,j), Ptmp, 1.0, Z(k,i));
02714           /* Z(j,j) is triangular: */
02715           Telement::trmm('L', 'U', false, -1.0, Z(j,j), Ptmp);
02716           Telement::add(1.0, Ptmp, Z(j,i));
02717         }
02718         Ptmp = A(i,i); /* Z(i,i) == I */
02719         for (int k = 0; k < i; k++) /* SYMMETRY USED */
02720           Telement::gemm_upper_tr_only(true, false, 1.0, 
02721                                        A(k,i), Z(k,i), 1.0, Ptmp);
02722         /* Z(i,i) == I !!                                                */
02723         /* Z(:,i) *= INCH(Ptmp) */
02724         Telement::syInch(Ptmp, Z(i,i), threshold, looking, version);
02725         for (int k = 0; k < i; k++) {         
02726           Telement::trmm('R', 'U', false, 1.0, Z(i,i), Z(k,i));
02727           /* INCH(Ptmp) is upper triangular*/
02728         } 
02729       }
02730     } /* end if left-looking inch */
02731     else { /* RIGHT-LOOKING INCH */
02732       assert(looking == right);  /* right and unstable */
02733       Telement Ptmp;
02734       Treal newThresh = 0;
02735       if (myThresh != 0 && Z.ncols() > 1)
02736         newThresh = myThresh * 0.0001;
02737       else
02738         newThresh = myThresh;
02739 
02740       for (int i = 0; i < Z.ncols(); i++) {
02741         /* Ptmp =  A(i,:) * Z(:,i)   */
02742         Ptmp = A(i,i); /* Z(i,i) == I */
02743         for (int k = 0; k < i; k++)
02744           Telement::gemm_upper_tr_only(true, false, 1.0, /* SYMMETRY USED */ 
02745                                        A(k,i), Z(k,i), 1.0, Ptmp);
02746         
02747         /* Z(:,i) *= INCH(Ptmp)  */
02748         Telement::syInch(Ptmp, Z(i,i), newThresh, looking, version);
02749         for (int k = 0; k < i; k++) {         
02750           Telement::trmm('R', 'U', false, 1.0, Z(i,i), Z(k,i));
02751           /* INCH(Ptmp) is upper triangular */
02752         }
02753           
02754         for (int j = i + 1; j < Z.ncols(); j++) {
02755           /* Compute Ptmp = Z(i,i)^T * A(i,:) * Z(:,j)  */
02756           /* Z(k,j) is nonzero for k = 0, 1, ...,i - 2, i - 1, j         */
02757           /* and Z(j,j) = I                                              */
02758           Ptmp = A(i,j);                  /* (Z(j,j) == I)             */
02759           for (int k = 0; k < i; k++)       /* Ptmp = A(i,:) * Z(:,j)  */
02760             Telement::gemm(true, false, 1.0, /* SYMMETRY USED */
02761                            A(k,i), Z(k,j), 1.0, Ptmp);
02762           Telement::trmm('L', 'U', true, 1.0, Z(i,i), Ptmp);
02763             
02764           for (int k = 0; k < i; k++)       /* Z(:,j) -= Z(:,i) * Ptmp   */
02765             Telement::gemm(false, false, -1.0,
02766                            Z(k,i), Ptmp, 1.0, Z(k,j));
02767           /* Z(i,i) is triangular: */
02768           Telement::trmm('L', 'U', false, -1.0, Z(i,i), Ptmp);
02769           Telement::add(1.0, Ptmp, Z(i,j));
02770         } /* end for j */
02771         
02772         /* Truncation starts here */
02773         if (threshold != 0) {
02774           for (int k = 0; k < i; k++) 
02775             Z(k,i).frob_thresh(myThresh);
02776         }
02777       } /* end for i */
02778     } /* end else right-looking inch */
02779   }
02780 
02781   template<class Treal, class Telement>
02782     void Matrix<Treal, Telement>:: 
02783     gersgorin(Treal& lmin, Treal& lmax) const {
02784     assert(!this->is_empty());
02785     if (this->nScalarsRows() == this->nScalarsCols()) {
02786       int size = this->nScalarsCols();
02787       Treal* colsums = new Treal[size];
02788       Treal* diag    = new Treal[size];
02789       for (int ind = 0; ind < size; ind++)
02790         colsums[ind] = 0;
02791       this->add_abs_col_sums(colsums);
02792       this->get_diagonal(diag);
02793       Treal tmp1 = colsums[0] - template_blas_fabs(diag[0]);
02794       Treal tmp2;
02795       lmin = diag[0] - tmp1;
02796       lmax = diag[0] + tmp1;
02797       for (int col = 1; col < size; col++) {
02798         tmp1 = colsums[col] - template_blas_fabs(diag[col]);
02799         tmp2 = diag[col] - tmp1;
02800         lmin = (tmp2 < lmin ? tmp2 : lmin);
02801         tmp2 = diag[col] + tmp1;
02802         lmax = (tmp2 > lmax ? tmp2 : lmax);
02803       }
02804       delete[] diag;
02805       delete[] colsums;
02806     }
02807     else
02808       throw Failure("Matrix<Treal, Telement, int>::gersgorin(Treal, Treal): "
02809                     "Matrix must be quadratic");
02810   }
02811 
02812 
02813   template<class Treal, class Telement> 
02814     void Matrix<Treal, Telement>::
02815     add_abs_col_sums(Treal* sums) const {
02816     assert(sums);
02817     if (!this->is_zero()) {
02818       int offset = 0;
02819       for (int col = 0; col < this->ncols(); col++) {   
02820         for (int row = 0; row < this->nrows(); row++) {
02821           (*this)(row,col).add_abs_col_sums(&sums[offset]);
02822         }
02823         offset += (*this)(0,col).nScalarsCols();
02824       }
02825     }  
02826   }
02827 
02828   template<class Treal, class Telement> 
02829     void Matrix<Treal, Telement>::
02830     get_diagonal(Treal* diag) const {
02831     assert(diag);
02832     assert(this->nrows() == this->ncols());
02833     if (this->is_zero()) 
02834       for (int rc = 0; rc < this->nScalarsCols(); rc++)
02835         diag[rc] = 0;
02836     else {
02837       int offset = 0;
02838       for (int rc = 0; rc < this->ncols(); rc++) {
02839         (*this)(rc,rc).get_diagonal(&diag[offset]);
02840         offset += (*this)(rc,rc).nScalarsCols();
02841       }
02842     }
02843   }
02844 
02845   template<class Treal, class Telement>
02846     size_t Matrix<Treal, Telement>::memory_usage() const {
02847     assert(!this->is_empty());
02848     size_t sum = 0; 
02849     if (this->is_zero())
02850       return (size_t)0;
02851     for (int ind = 0; ind < this->nElements(); ind++)
02852       sum += this->elements[ind].memory_usage();
02853     return sum;
02854   }
02855 
02856   template<class Treal, class Telement>
02857     size_t Matrix<Treal, Telement>::nnz() const {
02858     size_t sum = 0;
02859     if (!this->is_zero()) {
02860       for (int ind = 0; ind < this->nElements(); ind++)
02861         sum += this->elements[ind].nnz();
02862     }
02863     return sum;
02864   }
02865   template<class Treal, class Telement>
02866     size_t Matrix<Treal, Telement>::sy_nnz() const {
02867     size_t sum = 0;    
02868     if (!this->is_zero()) {
02869       /* Above diagonal */
02870       for (int col = 1; col < this->ncols(); col++)
02871         for (int row = 0; row < col; row++)
02872           sum += (*this)(row, col).nnz();
02873       /* Below diagonal */
02874       sum *= 2;
02875       /* Diagonal */
02876       for (int rc = 0; rc < this->nrows(); rc++)
02877         sum += (*this)(rc,rc).sy_nnz();
02878     }
02879     return sum;
02880   }
02881 
02882   template<class Treal, class Telement>
02883     size_t Matrix<Treal, Telement>::sy_nvalues() const {
02884     size_t sum = 0;
02885     if (!this->is_zero()) {
02886       /* Above diagonal */
02887       for (int col = 1; col < this->ncols(); col++)
02888         for (int row = 0; row < col; row++)
02889           sum += (*this)(row, col).nvalues();
02890       /* Diagonal */
02891       for (int rc = 0; rc < this->nrows(); rc++)
02892         sum += (*this)(rc,rc).sy_nvalues();
02893     }
02894     return sum;
02895   }
02896 
02897 
02898 
02899 
02900 
02901   /***************************************************************************/
02902   /***************************************************************************/
02903   /*           Specialization for Telement = Treal                           */
02904   /***************************************************************************/
02905   /***************************************************************************/
02906 
02907   template<class Treal>
02908     class Matrix<Treal>: public MatrixHierarchicBase<Treal> {
02909    
02910     public:
02911     typedef Vector<Treal, Treal> VectorType;
02912     friend class Vector<Treal, Treal>;
02913     
02914     Matrix()
02915       :MatrixHierarchicBase<Treal>() {
02916     }
02917       /*    Matrix(const Atomblock<Treal>& row_atoms, 
02918             const Atomblock<Treal>& col_atoms,
02919             bool z = true, int nr = 0, int nc = 0, char tr = 'N')
02920             :MatrixHierarchicBase<Treal>(row_atoms, col_atoms, z, nr, nc,tr) {}
02921       */
02922       
02923     void allocate() {
02924       assert(!this->is_empty());
02925       assert(this->is_zero());
02926       this->elements = new Treal[this->nElements()];
02927       for (int ind = 0; ind < this->nElements(); ++ind)
02928         this->elements[ind] = 0;
02929     }
02930     
02931     /* Full matrix assigns etc */
02932     void assignFromFull(std::vector<Treal> const & fullMat);
02933     void fullMatrix(std::vector<Treal> & fullMat) const; 
02934     void syFullMatrix(std::vector<Treal> & fullMat) const;
02935     void syUpTriFullMatrix(std::vector<Treal> & fullMat) const;
02936     
02937     /* Sparse matrix assigns etc */
02938     void assignFromSparse(std::vector<int> const & rowind, 
02939                           std::vector<int> const & colind, 
02940                           std::vector<Treal> const & values);
02941     void assignFromSparse(std::vector<int> const & rowind, 
02942                           std::vector<int> const & colind, 
02943                           std::vector<Treal> const & values,
02944                           std::vector<int> const & indexes);
02945     
02946     /* Adds values (+=) to elements */
02947     void addValues(std::vector<int> const & rowind, 
02948                    std::vector<int> const & colind, 
02949                    std::vector<Treal> const & values);
02950     void addValues(std::vector<int> const & rowind, 
02951                    std::vector<int> const & colind, 
02952                    std::vector<Treal> const & values,
02953                    std::vector<int> const & indexes);
02954     
02955     void syAssignFromSparse(std::vector<int> const & rowind, 
02956                             std::vector<int> const & colind, 
02957                             std::vector<Treal> const & values);
02958     
02959     void syAddValues(std::vector<int> const & rowind, 
02960                      std::vector<int> const & colind, 
02961                      std::vector<Treal> const & values);
02962     
02963     void getValues(std::vector<int> const & rowind, 
02964                    std::vector<int> const & colind, 
02965                    std::vector<Treal> & values) const;
02966     void getValues(std::vector<int> const & rowind, 
02967                    std::vector<int> const & colind, 
02968                    std::vector<Treal> & values,
02969                    std::vector<int> const & indexes) const;
02970     void syGetValues(std::vector<int> const & rowind, 
02971                      std::vector<int> const & colind, 
02972                      std::vector<Treal> & values) const;
02973     
02974     void getAllValues(std::vector<int> & rowind, 
02975                      std::vector<int> & colind, 
02976                      std::vector<Treal> & values) const;
02977     void syGetAllValues(std::vector<int> & rowind, 
02978                         std::vector<int> & colind, 
02979                         std::vector<Treal> & values) const;
02980 
02981     Matrix<Treal>& 
02982       operator=(const Matrix<Treal>& mat) {
02983       MatrixHierarchicBase<Treal>::operator=(mat);
02984       return *this;
02985     } 
02986 
02987     void clear(); 
02988     ~Matrix() {
02989       clear();
02990     }
02991 
02992     void writeToFile(std::ofstream & file) const;
02993     void readFromFile(std::ifstream & file);
02994 
02995     void random();
02996     void syRandom();
02997     void randomZeroStructure(Treal probabilityBeingZero);
02998     void syRandomZeroStructure(Treal probabilityBeingZero);
02999     template<typename TRule>
03000       void setElementsByRule(TRule & rule);
03001     template<typename TRule>
03002       void sySetElementsByRule(TRule & rule);
03003     
03004 
03005     void addIdentity(Treal alpha);        /* C += alpha * I */
03006 
03007     static void transpose(Matrix<Treal> const & A,
03008                           Matrix<Treal> & AT);
03009 
03010     void symToNosym();
03011     void nosymToSym();
03012     
03013     /* Set matrix to zero (k = 0) or identity (k = 1)                      */
03014     Matrix<Treal>& operator=(int const k);
03015 
03016     Matrix<Treal>& operator*=(const Treal alpha);
03017     
03018     static void gemm(const bool tA, const bool tB, const Treal alpha, 
03019                      const Matrix<Treal>& A, 
03020                      const Matrix<Treal>& B, 
03021                      const Treal beta, 
03022                      Matrix<Treal>& C);
03023     static void symm(const char side, const char uplo, const Treal alpha, 
03024                      const Matrix<Treal>& A, 
03025                      const Matrix<Treal>& B, 
03026                      const Treal beta, 
03027                      Matrix<Treal>& C);
03028     static void syrk(const char uplo, const bool tA, const Treal alpha, 
03029                      const Matrix<Treal>& A, 
03030                      const Treal beta, 
03031                      Matrix<Treal>& C);
03032     /* C = beta * C + alpha * A * A where A and C are symmetric */
03033     static void sysq(const char uplo, const Treal alpha, 
03034                      const Matrix<Treal>& A, 
03035                      const Treal beta, 
03036                      Matrix<Treal>& C);
03037     /* C = alpha * A * B + beta * C where A and B are symmetric */
03038     static void ssmm(const Treal alpha, 
03039                      const Matrix<Treal>& A, 
03040                      const Matrix<Treal>& B, 
03041                      const Treal beta, 
03042                      Matrix<Treal>& C);
03043     /* C = alpha * A * B + beta * C where A and B are symmetric
03044      * and only the upper triangle of C is computed.
03045      */
03046     static void ssmm_upper_tr_only(const Treal alpha, 
03047                                    const Matrix<Treal>& A, 
03048                                    const Matrix<Treal>& B, 
03049                                    const Treal beta, 
03050                                    Matrix<Treal>& C);
03051 
03052     static void trmm(const char side, const char uplo, const bool tA, 
03053                      const Treal alpha, 
03054                      const Matrix<Treal>& A, 
03055                      Matrix<Treal>& B);
03056 
03057     /* Frobenius norms */
03058     inline Treal frob() const {return template_blas_sqrt(frobSquared());} 
03059     Treal frobSquared() const;
03060     inline Treal syFrob() const {
03061       return template_blas_sqrt(this->syFrobSquared());
03062     }
03063     Treal syFrobSquared() const;
03064     
03065     inline static Treal frobDiff
03066       (const Matrix<Treal>& A,
03067        const Matrix<Treal>& B) {
03068       return template_blas_sqrt(frobSquaredDiff(A, B));
03069     }
03070     static Treal frobSquaredDiff
03071       (const Matrix<Treal>& A,
03072        const Matrix<Treal>& B);
03073     
03074     inline static Treal syFrobDiff
03075       (const Matrix<Treal>& A,
03076        const Matrix<Treal>& B) {
03077       return template_blas_sqrt(syFrobSquaredDiff(A, B));       
03078     }
03079     static Treal syFrobSquaredDiff
03080       (const Matrix<Treal>& A,
03081        const Matrix<Treal>& B);      
03082 
03083     Treal trace() const;
03084     static Treal trace_ab(const Matrix<Treal>& A,
03085                           const Matrix<Treal>& B); 
03086     static Treal trace_aTb(const Matrix<Treal>& A,
03087                            const Matrix<Treal>& B); 
03088     static Treal sy_trace_ab(const Matrix<Treal>& A,
03089                              const Matrix<Treal>& B); 
03090 
03091     static void add(const Treal alpha,   /* B += alpha * A */
03092                     const Matrix<Treal>& A, 
03093                     Matrix<Treal>& B);
03094     void assign(Treal const  alpha,   /* *this = alpha * A */
03095                 Matrix<Treal> const & A);
03096     
03097 
03098     /********** Help functions for thresholding */
03099     //    int getnnzBlocksLowestLevel() const;
03100     void getFrobSqLowestLevel(std::vector<Treal> & frobsq) const;
03101     void frobThreshLowestLevel
03102       (Treal const threshold, Matrix<Treal> * ErrorMatrix);
03103 
03104     void getFrobSqElementLevel(std::vector<Treal> & frobsq) const;
03105     void frobThreshElementLevel
03106       (Treal const threshold, Matrix<Treal> * ErrorMatrix);
03107     
03108 
03109 #if 0
03110     inline void frobThreshLowestLevel
03111       (Treal const threshold, 
03112        Matrix<Treal> * ErrorMatrix) {
03113       bool a,b;
03114       frobThreshLowestLevel(threshold, ErrorMatrix, a, b);
03115     }
03116 #endif
03117 
03118     void assignFrobNormsLowestLevel
03119       ( Matrix<Treal, Matrix<Treal> > const & A ) {
03120       if (!A.is_zero()) {
03121         if ( this->is_zero() )
03122           this->allocate();
03123         assert( this->nElements() == A.nElements() );
03124         for (int ind = 0; ind < this->nElements(); ind++)
03125           this->elements[ind] = A[ind].frob();
03126       }
03127       else
03128         this->clear();
03129     } 
03130 
03131     void syAssignFrobNormsLowestLevel( Matrix<Treal, Matrix<Treal> > const & A ) {
03132       if (!A.is_zero()) {
03133         if ( this->is_zero() )
03134           this->allocate();
03135         assert( this->nElements() == A.nElements() );
03136         for (int col = 1; col < this->ncols(); col++)
03137           for (int row = 0; row < col; row++)
03138             (*this)(row,col) = A(row,col).frob();
03139         for (int rc = 0; rc < this->nrows(); rc++)
03140           (*this)(rc,rc) = A(rc,rc).syFrob();
03141       }
03142       else
03143         this->clear();
03144     }
03145 
03146     void assignDiffFrobNormsLowestLevel( Matrix<Treal, Matrix<Treal> > const & A,
03147                                          Matrix<Treal, Matrix<Treal> > const & B ) {
03148       if ( A.is_zero() && B.is_zero() ) {
03149         // Both A and B are zero
03150         this->clear();
03151         return;
03152       }
03153       // At least one of A and B is nonzero
03154       if ( this->is_zero() )
03155         this->allocate();
03156       if ( !A.is_zero() && !B.is_zero() ) {
03157         // Both are nonzero
03158         assert( this->nElements() == A.nElements() );
03159         assert( this->nElements() == B.nElements() );
03160         for (int ind = 0; ind < this->nElements(); ind++) {
03161           Matrix<Treal> Diff(A[ind]);
03162           Matrix<Treal>::add( -1.0, B[ind], Diff );
03163           this->elements[ind] = Diff.frob();
03164         }
03165         return;
03166       }
03167       if ( !A.is_zero() ) {
03168         // A is nonzero
03169         this->assignFrobNormsLowestLevel( A );
03170         return;
03171       }
03172       if ( !B.is_zero() ) {
03173         // B is nonzero
03174         this->assignFrobNormsLowestLevel( B );
03175         return;
03176       }
03177     } 
03178     void syAssignDiffFrobNormsLowestLevel( Matrix<Treal, Matrix<Treal> > const & A,
03179                                            Matrix<Treal, Matrix<Treal> > const & B ) {
03180       if ( A.is_zero() && B.is_zero() ) {
03181         // Both A and B are zero
03182         this->clear();
03183         return;
03184       }
03185       // At least one of A and B is nonzero
03186       if ( this->is_zero() )
03187         this->allocate();
03188       if ( !A.is_zero() && !B.is_zero() ) {
03189         // Both are nonzero
03190         assert( this->nElements() == A.nElements() );
03191         assert( this->nElements() == B.nElements() );
03192         for (int col = 1; col < this->ncols(); col++)
03193           for (int row = 0; row < col; row++) {
03194             Matrix<Treal> Diff(A(row,col));
03195             Matrix<Treal>::add( -1.0, B(row,col), Diff );
03196             (*this)(row, col) = Diff.frob();
03197           }
03198         for (int rc = 0; rc < this->ncols(); rc++) {
03199           Matrix<Treal> Diff( A(rc,rc) );
03200           Matrix<Treal>::add( -1.0, B(rc,rc), Diff );
03201           (*this)(rc, rc) = Diff.syFrob();
03202         }
03203         return;
03204       }
03205       if ( !A.is_zero() ) {
03206         // A is nonzero
03207         this->syAssignFrobNormsLowestLevel( A );
03208         return;
03209       }
03210       if ( !B.is_zero() ) {
03211         // B is nonzero
03212         this->syAssignFrobNormsLowestLevel( B );
03213         return;
03214       }
03215     }
03216 
03217 
03218     void truncateAccordingToSparsityPattern( Matrix<Treal, Matrix<Treal> > & A ) const {
03219       if ( this->is_zero() ) 
03220         A.clear();
03221       else {
03222         assert( !A.is_zero() );
03223         assert( this->nElements() == A.nElements() );
03224         for (int ind = 0; ind < this->nElements(); ind++) 
03225           if (this->elements[ind] == 0)
03226             A[ind].clear();
03227       } 
03228     }
03229     
03230 
03231     /********** End of help functions for thresholding */
03232 
03233     static void gemm_upper_tr_only(const bool tA, const bool tB, 
03234                                    const Treal alpha, 
03235                                    const Matrix<Treal>& A, 
03236                                    const Matrix<Treal>& B, 
03237                                    const Treal beta, 
03238                                    Matrix<Treal>& C);
03239     static void sytr_upper_tr_only(char const side, const Treal alpha,
03240                                    Matrix<Treal>& A,
03241                                    const Matrix<Treal>& Z);
03242     static void trmm_upper_tr_only(const char side, const char uplo, 
03243                                    const bool tA, const Treal alpha, 
03244                                    const Matrix<Treal>& A, 
03245                                    Matrix<Treal>& B);
03246     static void trsytriplemm(char const side, 
03247                              const Matrix<Treal>& Z,
03248                              Matrix<Treal>& A);
03249 
03250     inline Treal frob_thresh(Treal const threshold,
03251                              Matrix<Treal> * ErrorMatrix = 0) {     
03252       return template_blas_sqrt
03253         (frob_squared_thresh(threshold * threshold, ErrorMatrix));
03254     }   
03255     /* Returns the Frobenius norm of the introduced error */
03256     
03257     Treal frob_squared_thresh(Treal const threshold,
03258                               Matrix<Treal> * ErrorMatrix = 0);
03259     
03260 
03261     static void inch(const Matrix<Treal>& A,
03262                      Matrix<Treal>& Z,
03263                      const Treal threshold = 0,
03264                      const side looking = left,
03265                      const inchversion version = unstable);
03266     static void syInch(const Matrix<Treal>& A,
03267                         Matrix<Treal>& Z,
03268                         const Treal threshold = 0, 
03269                         const side looking = left,
03270                         const inchversion version = unstable) {
03271       inch(A, Z, threshold, looking, version);
03272     }
03273 
03274     void gersgorin(Treal& lmin, Treal& lmax) const;
03275     void sy_gersgorin(Treal& lmin, Treal& lmax) const {
03276       Matrix<Treal> tmp = (*this);
03277       tmp.symToNosym();
03278       tmp.gersgorin(lmin, lmax);
03279       return;
03280     }
03281     void add_abs_col_sums(Treal* abscolsums) const;
03282     void get_diagonal(Treal* diag) const; /* Copy diagonal to the diag array */
03283     
03284     inline size_t memory_usage() const { /* Returns memory usage in bytes */
03285       assert(!this->is_empty());
03286       if (this->is_zero())
03287         return (size_t)0;
03288       else
03289         return (size_t)this->nElements() * sizeof(Treal);
03290     }
03291 
03292     inline size_t nnz() const {
03293       if (this->is_zero())
03294         return 0;
03295       else
03296         return this->nElements();
03297     } 
03298     inline size_t sy_nnz() const {
03299       if (this->is_zero())
03300         return 0;
03301       else
03302         return this->nElements();
03303     } 
03307     inline size_t nvalues() const {
03308       return nnz();
03309     } 
03312     size_t sy_nvalues() const {
03313       assert(this->nScalarsRows() == this->nScalarsCols());
03314       int n = this->nrows();
03315       return this->is_zero() ? 0 : n * (n + 1) / 2; 
03316     } 
03321     template<class Top> 
03322       Treal syAccumulateWith(Top & op) {
03323       Treal res = 0;
03324       if (!this->is_zero()) {
03325         int rowOffset = this->rows.getOffset();
03326         int colOffset = this->cols.getOffset();
03327         for (int col = 0; col < this->ncols(); col++) {
03328           for (int row = 0; row < col; row++) {
03329             res += 2*op.accumulate((*this)(row, col),
03330                                    rowOffset + row,
03331                                    colOffset + col);
03332           }
03333           res += op.accumulate((*this)(col, col),
03334                                rowOffset + col,
03335                                colOffset + col);
03336         }
03337       }
03338       return res;
03339     }
03340     template<class Top>
03341       Treal geAccumulateWith(Top & op) {
03342       Treal res = 0;
03343       if (!this->is_zero()) {
03344         int rowOffset = this->rows.getOffset();
03345         int colOffset = this->cols.getOffset();
03346         for (int col = 0; col < this->ncols(); col++)
03347           for (int row = 0; row < this->nrows(); row++)
03348             res += op.accumulate((*this)(row, col),
03349                                  rowOffset + row,
03350                                  colOffset + col);
03351       }
03352       return res;
03353     }
03354     
03355     static inline unsigned int level() {return 0;}
03356     
03357     Treal maxAbsValue() const {
03358       if (this->is_zero())
03359         return 0;
03360       else {
03361         Treal maxAbsGlobal = 0;
03362         Treal maxAbsLocal  = 0;
03363         for (int ind = 0; ind < this->nElements(); ++ind) {
03364           maxAbsLocal = template_blas_fabs(this->elements[ind]);
03365           maxAbsGlobal = maxAbsGlobal > maxAbsLocal ? 
03366             maxAbsGlobal : maxAbsLocal; 
03367         } /* end for */
03368         return maxAbsGlobal;
03369       } 
03370     }
03371     
03372     /* New routines above */    
03373     
03374 #if 0 /* OLD ROUTINES */   
03375 
03376 
03377 #if 0
03378     inline Matrix<Treal>& operator=(const Matrix<Treal>& mat) {
03379       this->MatrixHierarchicBase<Treal>::operator=(mat);
03380       std::cout<<"operator="<<std::endl;
03381     }
03382 #endif
03383 
03384 
03385 
03386 
03387     
03388 
03389 
03390 
03391 
03392 
03393     void trim_memory_usage();
03394 #if 1
03395     void write_to_buffer_count(int& zb_length, int& vb_length) const;
03396     void write_to_buffer(int* zerobuf, const int zb_length, 
03397                          Treal* valuebuf, const int vb_length,
03398                          int& zb_index, int& vb_index) const;
03399     void read_from_buffer(int* zerobuf, const int zb_length, 
03400                           Treal* valuebuf, const int vb_length,
03401                           int& zb_index, int& vb_index);
03402 #endif
03403 
03404 
03405 
03406 
03407     
03408     
03409     /* continue here */
03410 
03411     
03412 
03413     
03414 
03415 
03416 
03417     
03418     inline bool lowestLevel() const {return true;}
03419     //    inline unsigned int level() const {return 0;}
03420 
03421 #endif     /* END OF OLD ROUTINES */   
03422   protected:
03423   private:
03424     static const Treal ZERO;
03425     static const Treal ONE;
03426   }; /* end class specialization Matrix<Treal> */
03427 
03428   template<class Treal>
03429     const Treal Matrix<Treal>::ZERO = 0;
03430   template<class Treal>
03431     const Treal Matrix<Treal>::ONE = 1;
03432 
03433 #if 1
03434   /* Full matrix assigns etc */
03435   template<class Treal> 
03436     void Matrix<Treal>::
03437     assignFromFull(std::vector<Treal> const & fullMat) {
03438     int nTotalRows = this->rows.getNTotalScalars();
03439     int nTotalCols = this->cols.getNTotalScalars();
03440     assert((int)fullMat.size() == nTotalRows * nTotalCols);
03441     int rowOffset = this->rows.getOffset();
03442     int colOffset = this->cols.getOffset();
03443     if (this->is_zero())
03444       allocate();
03445     for (int col = 0; col < this->ncols(); col++) 
03446       for (int row = 0; row < this->nrows(); row++) 
03447         (*this)(row, col) = 
03448           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)];
03449   }
03450 
03451   template<class Treal> 
03452     void Matrix<Treal>::
03453     fullMatrix(std::vector<Treal> & fullMat) const {
03454     int nTotalRows = this->rows.getNTotalScalars();
03455     int nTotalCols = this->cols.getNTotalScalars();
03456     fullMat.resize(nTotalRows * nTotalCols);
03457     int rowOffset = this->rows.getOffset();
03458     int colOffset = this->cols.getOffset();
03459     if (this->is_zero()) {
03460       for (int col = 0; col < this->nScalarsCols(); col++)
03461         for (int row = 0; row < this->nScalarsRows(); row++) 
03462           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
03463     } 
03464     else {
03465       for (int col = 0; col < this->ncols(); col++) 
03466         for (int row = 0; row < this->nrows(); row++) 
03467           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 
03468             (*this)(row, col);
03469     }
03470   }
03471 
03472   template<class Treal> 
03473     void Matrix<Treal>::
03474     syFullMatrix(std::vector<Treal> & fullMat) const {
03475     int nTotalRows = this->rows.getNTotalScalars();
03476     int nTotalCols = this->cols.getNTotalScalars();
03477     fullMat.resize(nTotalRows * nTotalCols);
03478     int rowOffset = this->rows.getOffset();
03479     int colOffset = this->cols.getOffset();
03480     if (this->is_zero()) {
03481       for (int col = 0; col < this->nScalarsCols(); col++)
03482         for (int row = 0; row <= col; row++) {
03483           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
03484           fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
03485         }
03486     } 
03487     else {
03488       for (int col = 0; col < this->ncols(); col++) 
03489         for (int row = 0; row <= col; row++) {
03490           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 
03491             (*this)(row, col);
03492           fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 
03493             (*this)(row, col);
03494         }
03495     }
03496   }
03497   
03498   template<class Treal> 
03499     void Matrix<Treal>::
03500     syUpTriFullMatrix(std::vector<Treal> & fullMat) const {
03501     int nTotalRows = this->rows.getNTotalScalars();
03502     int nTotalCols = this->cols.getNTotalScalars();
03503     fullMat.resize(nTotalRows * nTotalCols);
03504     int rowOffset = this->rows.getOffset();
03505     int colOffset = this->cols.getOffset();
03506     if (this->is_zero()) {
03507       for (int col = 0; col < this->nScalarsCols(); col++)
03508         for (int row = 0; row <= this->nScalarsRows(); row++) {
03509           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 0;
03510           fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 0;
03511         }
03512     } 
03513     else {
03514       for (int col = 0; col < this->ncols(); col++) 
03515         for (int row = 0; row < this->nrows(); row++) {
03516           fullMat[(colOffset + col) * nTotalRows + (rowOffset + row)] = 
03517             (*this)(row, col);
03518           fullMat[(rowOffset + row) * nTotalRows + (colOffset + col)] = 
03519             (*this)(row, col);
03520         }
03521     }
03522   }
03523   
03524 #endif
03525   
03526   template<class Treal>
03527     void Matrix<Treal>::
03528     assignFromSparse(std::vector<int> const & rowind, 
03529                      std::vector<int> const & colind, 
03530                      std::vector<Treal> const & values) {
03531     assert(rowind.size() == colind.size() &&
03532            rowind.size() == values.size());
03533     std::vector<int> indexes(values.size());
03534     for (int ind = 0; ind < values.size(); ++ind) {
03535       indexes[ind] = ind;
03536     }
03537     assignFromSparse(rowind, colind, values, indexes);
03538   }
03539 
03540   template<class Treal>
03541     void Matrix<Treal>::
03542     assignFromSparse(std::vector<int> const & rowind, 
03543                      std::vector<int> const & colind, 
03544                      std::vector<Treal> const & values,
03545                      std::vector<int> const & indexes) {
03546     if (indexes.empty()) {
03547       this->clear();
03548       return;
03549     }
03550     if (this->is_zero())
03551       allocate();
03552     for (int ind = 0; ind < this->nElements(); ++ind)
03553       this->elements[ind] = 0;
03554     std::vector<int>::const_iterator it;
03555     for ( it = indexes.begin(); it < indexes.end(); it++ ) {
03556       (*this)(this->rows.whichBlock( rowind[*it] ),
03557               this->cols.whichBlock( colind[*it] )) = values[*it];
03558     }
03559   }
03560 
03561   
03562   template<class Treal>
03563     void Matrix<Treal>::
03564     addValues(std::vector<int> const & rowind, 
03565               std::vector<int> const & colind, 
03566               std::vector<Treal> const & values) {
03567     assert(rowind.size() == colind.size() &&
03568            rowind.size() == values.size());
03569     std::vector<int> indexes(values.size());
03570     for (int ind = 0; ind < values.size(); ++ind) {
03571       indexes[ind] = ind;
03572     }
03573     addValues(rowind, colind, values, indexes);
03574   }
03575   
03576   template<class Treal>
03577     void Matrix<Treal>::
03578     addValues(std::vector<int> const & rowind, 
03579               std::vector<int> const & colind, 
03580               std::vector<Treal> const & values,
03581               std::vector<int> const & indexes) {
03582     if (indexes.empty())
03583       return;
03584     if (this->is_zero())
03585       allocate();
03586     std::vector<int>::const_iterator it;
03587     for ( it = indexes.begin(); it < indexes.end(); it++ ) {
03588       (*this)(this->rows.whichBlock( rowind[*it] ),
03589               this->cols.whichBlock( colind[*it] )) += values[*it];
03590     }
03591   }
03592     
03593   template<class Treal>
03594     void Matrix<Treal>::
03595     syAssignFromSparse(std::vector<int> const & rowind, 
03596                        std::vector<int> const & colind, 
03597                        std::vector<Treal> const & values) {
03598     assert(rowind.size() == colind.size() &&
03599            rowind.size() == values.size());
03600     bool upperTriangleOnly = true;
03601     for (int ind = 0; ind < values.size(); ++ind) {
03602       upperTriangleOnly = 
03603         upperTriangleOnly && colind[ind] >= rowind[ind];
03604     }
03605     if (!upperTriangleOnly)
03606       throw Failure("Matrix<Treal>::"
03607                     "syAddValues(std::vector<int> const &, "
03608                     "std::vector<int> const &, "
03609                     "std::vector<Treal> const &, int const): "
03610                     "Only upper triangle can contain elements when assigning "
03611                     "symmetric or triangular matrix ");
03612     assignFromSparse(rowind, colind, values);
03613   }
03614     
03615   template<class Treal>
03616     void Matrix<Treal>::
03617     syAddValues(std::vector<int> const & rowind, 
03618                 std::vector<int> const & colind, 
03619                 std::vector<Treal> const & values) {
03620     assert(rowind.size() == colind.size() &&
03621            rowind.size() == values.size());
03622     bool upperTriangleOnly = true;
03623     for (int ind = 0; ind < values.size(); ++ind) {
03624       upperTriangleOnly = 
03625         upperTriangleOnly && colind[ind] >= rowind[ind];
03626     }
03627     if (!upperTriangleOnly)
03628       throw Failure("Matrix<Treal>::"
03629                     "syAddValues(std::vector<int> const &, "
03630                     "std::vector<int> const &, "
03631                     "std::vector<Treal> const &, int const): "
03632                     "Only upper triangle can contain elements when assigning "
03633                     "symmetric or triangular matrix ");
03634     addValues(rowind, colind, values);
03635   }
03636     
03637   template<class Treal>
03638     void Matrix<Treal>::
03639     getValues(std::vector<int> const & rowind, 
03640               std::vector<int> const & colind, 
03641               std::vector<Treal> & values) const {
03642     assert(rowind.size() == colind.size());
03643     values.resize(rowind.size(),0);
03644     std::vector<int> indexes(rowind.size());
03645     for (int ind = 0; ind < rowind.size(); ++ind) {
03646       indexes[ind] = ind;
03647     }
03648     getValues(rowind, colind, values, indexes);
03649   }
03650   
03651   template<class Treal>
03652     void Matrix<Treal>::
03653     getValues(std::vector<int> const & rowind, 
03654               std::vector<int> const & colind, 
03655               std::vector<Treal> & values,
03656               std::vector<int> const & indexes) const {
03657     assert(!this->is_empty());
03658     if (indexes.empty())
03659       return;
03660     std::vector<int>::const_iterator it;
03661     if (this->is_zero()) {
03662       for ( it = indexes.begin(); it < indexes.end(); it++ )
03663         values[*it] = 0;
03664       return;
03665     }
03666     for ( it = indexes.begin(); it < indexes.end(); it++ ) 
03667       values[*it] = (*this)(this->rows.whichBlock( rowind[*it] ),
03668                             this->cols.whichBlock( colind[*it] ));
03669   }
03670   
03671 
03672   template<class Treal>
03673     void Matrix<Treal>::
03674     syGetValues(std::vector<int> const & rowind, 
03675                 std::vector<int> const & colind, 
03676                 std::vector<Treal> & values) const {
03677     assert(rowind.size() == colind.size());
03678     bool upperTriangleOnly = true;
03679     for (int ind = 0; ind < rowind.size(); ++ind) {
03680       upperTriangleOnly = 
03681         upperTriangleOnly && colind[ind] >= rowind[ind];
03682     }
03683     if (!upperTriangleOnly)
03684       throw Failure("Matrix<Treal>::"
03685                     "syGetValues(std::vector<int> const &, "
03686                     "std::vector<int> const &, "
03687                     "std::vector<Treal> const &, int const): "
03688                     "Only upper triangle when retrieving elements from "
03689                     "symmetric or triangular matrix ");
03690     getValues(rowind, colind, values);
03691   }
03692 
03693   template<class Treal>
03694     void Matrix<Treal>::
03695     getAllValues(std::vector<int> & rowind, 
03696                  std::vector<int> & colind, 
03697                  std::vector<Treal> & values) const {
03698     assert(rowind.size() == colind.size() &&
03699            rowind.size() == values.size());
03700     if (!this->is_zero()) 
03701       for (int col = 0; col < this->ncols(); col++)
03702         for (int row = 0; row < this->nrows(); row++) {
03703           rowind.push_back(this->rows.getOffset() + row);
03704           colind.push_back(this->cols.getOffset() + col);
03705           values.push_back((*this)(row, col));
03706         }
03707   }
03708   
03709   template<class Treal>
03710     void Matrix<Treal>::
03711     syGetAllValues(std::vector<int> & rowind, 
03712                    std::vector<int> & colind, 
03713                    std::vector<Treal> & values) const {
03714     assert(rowind.size() == colind.size() &&
03715            rowind.size() == values.size());
03716     if (!this->is_zero()) 
03717       for (int col = 0; col < this->ncols(); ++col) 
03718         for (int row = 0; row <= col; ++row) {
03719           rowind.push_back(this->rows.getOffset() + row);
03720           colind.push_back(this->cols.getOffset() + col);
03721           values.push_back((*this)(row, col));
03722         }
03723   }
03724 
03725     
03726   template<class Treal>
03727     void Matrix<Treal>::clear() {
03728     delete[] this->elements;
03729     this->elements = 0;
03730   }
03731 
03732   template<class Treal>
03733     void Matrix<Treal>:: 
03734     writeToFile(std::ofstream & file) const {
03735     int const ZERO = 0;
03736     int const ONE  = 1;
03737     if (this->is_zero()) {
03738       char * tmp = (char*)&ZERO;
03739       file.write(tmp,sizeof(int));
03740     }
03741     else {
03742       char * tmp = (char*)&ONE;
03743       file.write(tmp,sizeof(int));
03744       char * tmpel = (char*)this->elements;
03745       file.write(tmpel,sizeof(Treal) * this->nElements());
03746     }
03747   }
03748 
03749   template<class Treal>
03750     void Matrix<Treal>:: 
03751     readFromFile(std::ifstream & file) {
03752     int const ZERO = 0;
03753     int const ONE  = 1;
03754     char tmp[sizeof(int)];
03755     file.read(tmp, (std::ifstream::pos_type)sizeof(int));
03756     switch ((int)*tmp) {
03757     case ZERO:
03758       this->clear();
03759       break;
03760     case ONE:
03761       if (this->is_zero())
03762         allocate();
03763       file.read((char*)this->elements, sizeof(Treal) * this->nElements());
03764       break;
03765     default:
03766       throw Failure("Matrix<Treal>::" 
03767                     "readFromFile(std::ifstream & file):"
03768                     "File corruption, int value not 0 or 1");
03769     }
03770   }
03771 
03772   template<class Treal> 
03773     void Matrix<Treal>::random() {
03774     if (this->is_zero())
03775       allocate();
03776     for (int ind = 0; ind < this->nElements(); ind++)
03777       this->elements[ind] = rand() / (Treal)RAND_MAX;
03778   }
03779   
03780   template<class Treal> 
03781     void Matrix<Treal>::syRandom() {
03782     if (this->is_zero())
03783       allocate();
03784     /* Above diagonal */
03785     for (int col = 1; col < this->ncols(); col++)
03786       for (int row = 0; row < col; row++)
03787         (*this)(row, col) = rand() / (Treal)RAND_MAX;;
03788     /* Diagonal */
03789     for (int rc = 0; rc < this->nrows(); rc++)
03790       (*this)(rc,rc) = rand() / (Treal)RAND_MAX;;
03791   }
03792   
03793   template<class Treal> 
03794     void Matrix<Treal>::
03795     randomZeroStructure(Treal probabilityBeingZero) {
03796     if (!this->highestLevel() && 
03797         probabilityBeingZero > rand() / (Treal)RAND_MAX)
03798       this->clear();
03799     else 
03800       this->random();
03801   }
03802 
03803   template<class Treal> 
03804     void  Matrix<Treal>::
03805     syRandomZeroStructure(Treal probabilityBeingZero) {
03806     if (!this->highestLevel() && 
03807         probabilityBeingZero > rand() / (Treal)RAND_MAX)
03808       this->clear();
03809     else 
03810       this->syRandom();
03811   }
03812 
03813   template<class Treal> 
03814     template<typename TRule>
03815     void Matrix<Treal>::
03816     setElementsByRule(TRule & rule) {
03817     if (this->is_zero()) 
03818       allocate();
03819     for (int col = 0; col < this->ncols(); col++)
03820       for (int row = 0; row < this->nrows(); row++)
03821         (*this)(row,col) = rule.set(this->rows.getOffset() + row,
03822                                       this->cols.getOffset() + col);
03823   }
03824   template<class Treal> 
03825     template<typename TRule>
03826     void Matrix<Treal>::
03827     sySetElementsByRule(TRule & rule) {
03828     if (this->is_zero()) 
03829       allocate(); 
03830     /* Upper triangle */
03831     for (int col = 0; col < this->ncols(); col++)
03832       for (int row = 0; row <= col; row++)
03833         (*this)(row,col) = rule.set(this->rows.getOffset() + row,
03834                                       this->cols.getOffset() + col);
03835   }
03836 
03837 
03838   template<class Treal> 
03839     void Matrix<Treal>::
03840     addIdentity(Treal alpha) {
03841     if (this->is_empty())
03842       throw Failure("Matrix<Treal>::addIdentity(Treal): "
03843                     "Cannot add identity to empty matrix.");
03844     if (this->ncols() != this->nrows()) 
03845       throw Failure("Matrix<Treal, Telement>::addIdentity(Treal): "
03846                     "Matrix must be square to add identity");
03847     if (this->is_zero()) 
03848       allocate();
03849     for (int ind = 0; ind < this->ncols(); ind++)
03850       (*this)(ind,ind) += alpha;
03851   }
03852 
03853   template<class Treal> 
03854     void Matrix<Treal>::
03855     transpose(Matrix<Treal> const & A, Matrix<Treal> & AT) {
03856     if (A.is_zero()) { /* Condition also matches empty matrices. */
03857       AT.rows = A.cols;
03858       AT.cols = A.rows;
03859       delete[] AT.elements;
03860       AT.elements = 0;
03861       return;
03862     }
03863     if (AT.is_zero() || (AT.nElements() != A.nElements())) {
03864       delete[] AT.elements;
03865       AT.elements = new Treal[A.nElements()];
03866     }
03867     AT.cols = A.rows;
03868     AT.rows = A.cols;
03869     for (int row = 0; row < AT.nrows(); row++)
03870       for (int col = 0; col < AT.ncols(); col++)
03871         AT(row,col) = A(col,row);
03872   }
03873 
03874   template<class Treal> 
03875     void Matrix<Treal>::
03876     symToNosym() {
03877     if (this->nrows() == this->ncols()) {
03878       if (!this->is_zero()) {
03879         /* Diagonal should be fine */
03880         /* Fix the lower triangle */
03881         for (int row = 1; row < this->nrows(); row++)
03882           for (int col = 0; col < row; col++)
03883             (*this)(row, col) = (*this)(col, row);
03884       }
03885     }
03886     else
03887       throw Failure("Matrix<Treal>::symToNosym(): "
03888                     "Only quadratic matrices can be symmetric");
03889   }
03890 
03891   template<class Treal> 
03892     void Matrix<Treal>::
03893     nosymToSym() {
03894     if (this->nrows() == this->ncols()) {
03895       if (!this->is_zero()) {
03896         /* Diagonal should be fine */
03897         /* Fix the lower triangle */
03898         for (int col = 0; col < this->ncols() - 1; col++)
03899           for (int row = col + 1; row < this->nrows(); row++)
03900             (*this)(row, col) = 0;
03901       }
03902     }
03903     else
03904       throw Failure("Matrix<Treal>::nosymToSym(): "
03905                     "Only quadratic matrices can be symmetric");
03906   }
03907 
03908   template<class Treal>
03909     Matrix<Treal>& 
03910     Matrix<Treal>::operator=(int const k) {
03911     switch (k) {
03912     case 0:
03913       this->clear();
03914       break;
03915     case 1:
03916       if (this->ncols() != this->nrows())
03917         throw Failure("Matrix<Treal>::operator=(int k = 1): "
03918                       "Matrix must be quadratic to "
03919                       "become identity matrix.");
03920       this->clear();
03921       this->allocate();
03922       for (int rc = 0; rc < this->ncols(); rc++) /*Set diagonal to identity*/
03923         (*this)(rc,rc) = 1;
03924       break;
03925     default:
03926       throw Failure
03927         ("Matrix<Treal>::operator=(int k) only implemented for k = 0, k = 1");
03928     }
03929     return *this;
03930   }
03931 
03932   template<class Treal>
03933     Matrix<Treal>& Matrix<Treal>:: 
03934     operator*=(const Treal alpha) {
03935     if (!this->is_zero() && alpha != 1) {
03936       for (int ind = 0; ind < this->nElements(); ind++)
03937         this->elements[ind] *= alpha;
03938     }
03939     return *this;
03940   }
03941 
03942   template<class Treal> 
03943     void Matrix<Treal>::
03944     gemm(const bool tA, const bool tB, const Treal alpha, 
03945          const Matrix<Treal>& A, 
03946          const Matrix<Treal>& B, 
03947          const Treal beta, 
03948          Matrix<Treal>& C) {
03949     assert(!A.is_empty());
03950     assert(!B.is_empty());
03951     if (C.is_empty()) {
03952       assert(beta == 0);
03953       if (tA)
03954         C.resetRows(A.cols);
03955       else
03956         C.resetRows(A.rows);
03957       if (tB)
03958         C.resetCols(B.rows);
03959       else
03960         C.resetCols(B.cols);
03961     }      
03962     
03963     if ((A.is_zero() || B.is_zero() || alpha == 0) && 
03964         (C.is_zero() || beta == 0))
03965       C = 0;
03966     else {
03967       Treal beta_tmp = beta;
03968       if (C.is_zero()) {
03969         C.allocate();
03970         beta_tmp = 0;
03971       }
03972         
03973       if (!A.is_zero() && !B.is_zero() && alpha != 0) {
03974         if (!tA && !tB) 
03975           if (A.ncols() == B.nrows() && 
03976               A.nrows() == C.nrows() &&
03977               B.ncols() == C.ncols()) 
03978             mat::gemm("N", "N", &A.nrows(), &B.ncols(), &A.ncols(), &alpha,
03979                       A.elements, &A.nrows(), B.elements, &B.nrows(),
03980                       &beta_tmp, C.elements, &C.nrows());
03981           else 
03982             throw Failure("Matrix<Treal>::"
03983                           "gemm(bool, bool, Treal, Matrix<Treal>, "
03984                           "Matrix<Treal>, Treal, "
03985                           "Matrix<Treal>): "
03986                           "Incorrect matrixdimensions for multiplication");
03987         else if (tA && !tB)
03988           if (A.nrows() == B.nrows() &&
03989               A.ncols() == C.nrows() &&
03990               B.ncols() == C.ncols()) 
03991             mat::gemm("T", "N", &A.ncols(), &B.ncols(), &A.nrows(), &alpha,
03992                       A.elements, &A.nrows(), B.elements, &B.nrows(),
03993                       &beta_tmp, C.elements, &C.nrows());
03994           else
03995             throw Failure("Matrix<Treal>::"
03996                           "gemm(bool, bool, Treal, Matrix<Treal>, "
03997                           "Matrix<Treal>, Treal, "
03998                           "Matrix<Treal>): "
03999                           "Incorrect matrixdimensions for multiplication");
04000         else if (!tA && tB)
04001           if (A.ncols() == B.ncols() && 
04002               A.nrows() == C.nrows() &&
04003               B.nrows() == C.ncols()) 
04004             mat::gemm("N", "T", &A.nrows(), &B.nrows(), &A.ncols(), &alpha,
04005                       A.elements, &A.nrows(), B.elements, &B.nrows(),
04006                       &beta_tmp, C.elements, &C.nrows());
04007           else 
04008             throw Failure("Matrix<Treal>::"
04009                           "gemm(bool, bool, Treal, Matrix<Treal>, "
04010                           "Matrix<Treal>, Treal, "
04011                           "Matrix<Treal>): "
04012                           "Incorrect matrixdimensions for multiplication");
04013         else if (tA && tB)
04014           if (A.nrows() == B.ncols() && 
04015               A.ncols() == C.nrows() &&
04016               B.nrows() == C.ncols()) 
04017             mat::gemm("T", "T", &A.ncols(), &B.nrows(), &A.nrows(), &alpha,
04018                       A.elements, &A.nrows(), B.elements, &B.nrows(),
04019                       &beta_tmp, C.elements, &C.nrows());
04020           else 
04021             throw Failure("Matrix<Treal>::"
04022                           "gemm(bool, bool, Treal, Matrix<Treal>, "
04023                           "Matrix<Treal>, Treal, "
04024                           "Matrix<Treal>): "
04025                           "Incorrect matrixdimensions for multiplication");
04026         else throw Failure("Matrix<Treal>::"
04027                            "gemm(bool, bool, Treal, Matrix<Treal>, "
04028                            "Matrix<Treal>, Treal, "
04029                            "Matrix<Treal>):Very strange error!!");
04030       }
04031       else 
04032         C *= beta;
04033     }
04034   }
04035 
04036 
04037   template<class Treal> 
04038     void Matrix<Treal>::
04039     symm(const char side, const char uplo, const Treal alpha, 
04040          const Matrix<Treal>& A, 
04041          const Matrix<Treal>& B, 
04042          const Treal beta, 
04043          Matrix<Treal>& C) {
04044     assert(!A.is_empty());
04045     assert(!B.is_empty());
04046     assert(A.nrows() == A.ncols());
04047     //    int dimA = A.nrows();
04048     if (C.is_empty()) {
04049       assert(beta == 0);
04050       if (side =='L') {
04051         C.resetRows(A.rows);
04052         C.resetCols(B.cols);
04053       } 
04054       else {
04055         assert(side == 'R');
04056         C.resetRows(B.rows);
04057         C.resetCols(A.cols);
04058       }
04059     }      
04060     
04061     if ((A.is_zero() || B.is_zero() || alpha == 0) && 
04062         (C.is_zero() || beta == 0))
04063       C = 0;
04064     else {
04065       Treal beta_tmp = beta;
04066       if (C.is_zero()) {
04067         C.allocate();
04068         beta_tmp = 0;
04069       }
04070       if (!A.is_zero() && !B.is_zero() && alpha != 0) {
04071         if (A.nrows() == A.ncols() && C.nrows() == B.nrows() && 
04072             C.ncols() == B.ncols() && 
04073             ((side == 'L' && A.ncols() == B.nrows()) ||
04074              (side == 'R' && B.ncols() == A.nrows())))
04075           mat::symm(&side, &uplo, &C.nrows(), &C.ncols(), &alpha, 
04076                     A.elements, &A.nrows(), B.elements, &B.nrows(),
04077                     &beta_tmp, C.elements, &C.nrows());
04078         else
04079           throw Failure("Matrix<Treal>::symm: Incorrect matrix "
04080                         "dimensions for symmetric multiply");
04081       } /* end if (!A.is_zero() && !B.is_zero() && alpha != 0) */
04082       else 
04083         C *= beta;
04084     }
04085   }
04086 
04087   template<class Treal> 
04088     void Matrix<Treal>::
04089     syrk(const char uplo, const bool tA, const Treal alpha, 
04090          const Matrix<Treal>& A, 
04091          const Treal beta, 
04092          Matrix<Treal>& C) {
04093     assert(!A.is_empty());
04094     if (C.is_empty()) {
04095       assert(beta == 0);
04096       if (tA) {
04097         C.resetRows(A.cols);
04098         C.resetCols(A.cols);
04099       } 
04100       else {
04101         C.resetRows(A.rows);
04102         C.resetCols(A.rows);
04103       }
04104     }      
04105     if (C.nrows() == C.ncols() && 
04106         ((!tA && A.nrows() == C.nrows()) || (tA && A.ncols() == C.nrows())))
04107       if (alpha != 0 && !A.is_zero()) {
04108         Treal beta_tmp = beta;
04109         if (C.is_zero()) {
04110           C.allocate();
04111           beta_tmp = 0;
04112         }       
04113         if (!tA) {
04114           mat::syrk(&uplo, "N", &C.nrows(), &A.ncols(), &alpha, 
04115                     A.elements, &A.nrows(), &beta_tmp, 
04116                     C.elements, &C.nrows());
04117         } /* end if (!tA) */
04118         else
04119           mat::syrk(&uplo, "T", &C.nrows(), &A.nrows(), &alpha, 
04120                     A.elements, &A.nrows(), &beta_tmp, 
04121                     C.elements, &C.nrows());
04122       } 
04123       else
04124         C *= beta;
04125     else
04126       throw Failure("Matrix<Treal>::syrk: Incorrect matrix "
04127                     "dimensions for symmetric rank-k update");
04128   }
04129   
04130   template<class Treal> 
04131     void Matrix<Treal>::
04132     sysq(const char uplo, const Treal alpha, 
04133          const Matrix<Treal>& A, 
04134          const Treal beta, 
04135          Matrix<Treal>& C) {
04136     assert(!A.is_empty());
04137     if (C.is_empty()) {
04138       assert(beta == 0);
04139       C.resetRows(A.rows);
04140       C.resetCols(A.cols);
04141     }
04142     /* FIXME: slow copy */
04143     Matrix<Treal> tmpA = A; 
04144     tmpA.symToNosym();
04145     Matrix<Treal>::syrk(uplo, false, alpha, tmpA, beta, C);
04146   }
04147 
04148   /* C = alpha * A * B + beta * C where A and B are symmetric */
04149   template<class Treal> 
04150     void Matrix<Treal>::
04151     ssmm(const Treal alpha, 
04152          const Matrix<Treal>& A, 
04153          const Matrix<Treal>& B, 
04154          const Treal beta, 
04155          Matrix<Treal>& C) {
04156     assert(!A.is_empty());
04157     assert(!B.is_empty());
04158     if (C.is_empty()) {
04159       assert(beta == 0);
04160       C.resetRows(A.rows);
04161       C.resetCols(B.cols);
04162     }
04163     /* FIXME: slow copy */
04164     Matrix<Treal> tmpB = B; 
04165     tmpB.symToNosym();
04166     Matrix<Treal>::symm('L', 'U', alpha, A, tmpB, beta, C);
04167   }
04168   
04169   /* C = alpha * A * B + beta * C where A and B are symmetric
04170    * and only the upper triangle of C is computed.
04171    */
04172   template<class Treal> 
04173     void Matrix<Treal>::
04174     ssmm_upper_tr_only(const Treal alpha, 
04175                        const Matrix<Treal>& A, 
04176                        const Matrix<Treal>& B, 
04177                        const Treal beta, 
04178                        Matrix<Treal>& C) {
04179     /* FIXME: Symmetry is not utilized. */
04180     ssmm(alpha, A, B, beta, C);
04181     C.nosymToSym();   
04182   }
04183 
04184 
04185   template<class Treal> 
04186     void Matrix<Treal>::
04187     trmm(const char side, const char uplo, const bool tA, 
04188          const Treal alpha, 
04189          const Matrix<Treal>& A, 
04190          Matrix<Treal>& B) {
04191     assert(!A.is_empty());
04192     assert(!B.is_empty());
04193     if (alpha != 0 && !A.is_zero() && !B.is_zero())
04194       if (((side == 'R' && B.ncols() == A.nrows()) || 
04195            (side == 'L' && A.ncols() == B.nrows())) &&
04196           A.nrows() == A.ncols())
04197         if (!tA)
04198           mat::trmm(&side, &uplo, "N", "N", &B.nrows(), &B.ncols(),
04199                     &alpha, A.elements, &A.nrows(), B.elements, &B.nrows());
04200         else
04201           mat::trmm(&side, &uplo, "T", "N", &B.nrows(), &B.ncols(),
04202                     &alpha, A.elements, &A.nrows(), B.elements, &B.nrows());
04203       else
04204         throw Failure("Matrix<Treal>::trmm: "
04205                       "Incorrect matrix dimensions for multiplication");
04206     else
04207       B = 0;
04208   }
04209 
04210   template<class Treal>
04211     Treal Matrix<Treal>::frobSquared() const {
04212     assert(!this->is_empty());
04213     if (this->is_zero()) 
04214       return 0;
04215     else {
04216       Treal sum(0);
04217       for (int i = 0; i < this->nElements(); i++)
04218         sum += this->elements[i] * this->elements[i];
04219       return sum;
04220     }
04221   }
04222 
04223   template<class Treal>
04224     Treal Matrix<Treal>::
04225     syFrobSquared() const {
04226     assert(!this->is_empty());
04227     if (this->nrows() != this->ncols()) 
04228       throw Failure("Matrix<Treal>::syFrobSquared(): "
04229                     "Matrix must be have equal number of rows "
04230                     "and cols to be symmetric");
04231     Treal sum(0);
04232     if (!this->is_zero()) {
04233       for (int col = 1; col < this->ncols(); col++)
04234         for (int row = 0; row < col; row++)
04235           sum += 2 * (*this)(row, col) * (*this)(row, col);
04236       for (int rc = 0; rc < this->ncols(); rc++)
04237         sum += (*this)(rc, rc) * (*this)(rc, rc);
04238     }
04239     return sum;
04240   }
04241 
04242   template<class Treal>
04243     Treal Matrix<Treal>::
04244     frobSquaredDiff
04245     (const Matrix<Treal>& A,
04246      const Matrix<Treal>& B) {
04247     assert(!A.is_empty());
04248     assert(!B.is_empty());
04249     if (A.nrows() != B.nrows() || A.ncols() != B.ncols()) 
04250       throw Failure("Matrix<Treal>::frobSquaredDiff: "
04251                     "Incorrect matrix dimensions");
04252     Treal sum(0);
04253     if (!A.is_zero() && !B.is_zero()) {
04254       Treal diff;
04255       for (int i = 0; i < A.nElements(); i++) {
04256         diff = A.elements[i] - B.elements[i];
04257         sum += diff * diff;
04258       }
04259     }
04260     else if (A.is_zero() && !B.is_zero()) 
04261       sum = B.frobSquared();
04262     else if (!A.is_zero() && B.is_zero())
04263       sum = A.frobSquared();
04264     /* sum is already zero if A.is_zero() && B.is_zero() */
04265     return sum;
04266   }
04267 
04268   
04269   template<class Treal>
04270     Treal Matrix<Treal>::
04271     syFrobSquaredDiff
04272     (const Matrix<Treal>& A,
04273      const Matrix<Treal>& B) {
04274     assert(!A.is_empty());
04275     assert(!B.is_empty());
04276     if (A.nrows() != B.nrows() || 
04277         A.ncols() != B.ncols() || 
04278         A.nrows() != A.ncols()) 
04279       throw Failure("Matrix<Treal>::syFrobSquaredDiff: "
04280                     "Incorrect matrix dimensions");
04281     Treal sum(0);
04282     if (!A.is_zero() && !B.is_zero()) {
04283       Treal diff;
04284       for (int col = 1; col < A.ncols(); col++)
04285         for (int row = 0; row < col; row++) {
04286           diff = A(row, col) - B(row, col);
04287           sum += 2 * diff * diff;
04288         }
04289       for (int rc = 0; rc < A.ncols(); rc++) {
04290         diff = A(rc, rc) - B(rc, rc);
04291         sum += diff * diff;
04292       }
04293     }
04294     else if (A.is_zero() && !B.is_zero()) 
04295       sum = B.syFrobSquared();
04296     else if (!A.is_zero() && B.is_zero())
04297       sum = A.syFrobSquared();
04298     /* sum is already zero if A.is_zero() && B.is_zero() */
04299     return sum;
04300   }
04301 
04302   template<class Treal> 
04303     Treal Matrix<Treal>::
04304     trace() const {
04305     assert(!this->is_empty());
04306     if (this->nrows() != this->ncols())
04307       throw Failure("Matrix<Treal>::trace(): Matrix must be quadratic");  
04308     if (this->is_zero())
04309       return 0;
04310     else {
04311       Treal sum = 0;
04312       for (int rc = 0; rc < this->ncols(); rc++)
04313         sum += (*this)(rc,rc);
04314       return sum;
04315     }
04316   }
04317 
04318   template<class Treal> 
04319     Treal Matrix<Treal>::
04320     trace_ab(const Matrix<Treal>& A,
04321              const Matrix<Treal>& B) {
04322     assert(!A.is_empty());
04323     assert(!B.is_empty());
04324     if (A.nrows() != B.ncols() || A.ncols() != B.nrows()) 
04325       throw Failure("Matrix<Treal>::trace_ab: "
04326                     "Wrong matrix dimensions for trace(A * B)"); 
04327     Treal tr = 0;      
04328     if (!A.is_zero() && !B.is_zero())
04329       for (int rc = 0; rc < A.nrows(); rc++)
04330         for (int colA = 0; colA < A.ncols(); colA++)
04331           tr += A(rc,colA) * B(colA, rc);
04332     return tr;
04333   } 
04334   
04335   template<class Treal> 
04336     Treal  Matrix<Treal>::
04337     trace_aTb(const Matrix<Treal>& A,
04338               const Matrix<Treal>& B) {
04339     assert(!A.is_empty());
04340     assert(!B.is_empty());
04341     if (A.ncols() != B.ncols() || A.nrows() != B.nrows())
04342       throw Failure("Matrix<Treal>::trace_aTb: "
04343                     "Wrong matrix dimensions for trace(A' * B)"); 
04344     Treal tr = 0;
04345     if (!A.is_zero() && !B.is_zero())
04346       for (int rc = 0; rc < A.ncols(); rc++)
04347         for (int rowB = 0; rowB < B.nrows(); rowB++)
04348           tr += A(rowB,rc) * B(rowB, rc);
04349     return tr;
04350   }
04351 
04352   template<class Treal> 
04353     Treal Matrix<Treal>::
04354     sy_trace_ab(const Matrix<Treal>& A,
04355                 const Matrix<Treal>& B) {
04356     assert(!A.is_empty());
04357     assert(!B.is_empty());
04358     if (A.nrows() != B.ncols() || A.ncols() != B.nrows() || 
04359         A.nrows() != A.ncols()) 
04360       throw Failure("Matrix<Treal>::sy_trace_ab: "
04361                     "Wrong matrix dimensions for symmetric trace(A * B)");
04362     if (A.is_zero() || B.is_zero())
04363       return 0;
04364     /* Now we know both A and B are nonzero */
04365     Treal tr = 0;
04366     /* Diagonal first */
04367     for (int rc = 0; rc < A.nrows(); rc++)
04368       tr += A(rc,rc) * B(rc, rc);
04369     /* Using that trace of transpose is equal to that without transpose: */
04370     for (int colA = 1; colA < A.ncols(); colA++)
04371       for (int rowA = 0; rowA < colA; rowA++)
04372         tr += 2 * A(rowA, colA) * B(rowA, colA);
04373     return tr;
04374   }
04375 
04376   template<class Treal>
04377     void Matrix<Treal>:: 
04378     add(const Treal alpha,   /* B += alpha * A */
04379         const Matrix<Treal>& A, 
04380         Matrix<Treal>& B) {
04381     assert(!A.is_empty());
04382     assert(!B.is_empty());
04383     if (A.nrows() != B.nrows() || A.ncols() != B.ncols())
04384       throw Failure("Matrix<Treal>::add: "
04385                     "Wrong matrix dimensions for addition");
04386     if (!A.is_zero() && alpha != 0) {
04387       if (B.is_zero())
04388         B.allocate();
04389       for (int ind = 0; ind < A.nElements(); ind++)
04390         B.elements[ind] += alpha * A.elements[ind];
04391     }
04392   }
04393 
04394   template<class Treal>
04395     void Matrix<Treal>:: 
04396     assign(Treal const  alpha,   /* *this = alpha * A */
04397            Matrix<Treal> const & A) {
04398     assert(!A.is_empty());
04399     if (this->is_empty()) {
04400       this->resetRows(A.rows);
04401       this->resetCols(A.cols);
04402     }
04403     Matrix<Treal>::add(alpha, A, *this);
04404   }
04405   
04406 
04407   /********** Help functions for thresholding */
04408   
04409   template<class Treal>
04410     void Matrix<Treal>::
04411     getFrobSqLowestLevel(std::vector<Treal> & frobsq) const {
04412     if (!this->is_zero()) 
04413       frobsq.push_back(this->frobSquared());
04414   }
04415 
04416   template<class Treal>
04417     void Matrix<Treal>::
04418     getFrobSqElementLevel(std::vector<Treal> & frobsq) const {
04419     if (!this->is_zero()) 
04420       for (int ind = 0; ind < this->nElements(); ind++) 
04421         if ( this->elements[ind] != 0 ) // Add nonzero elements only
04422           frobsq.push_back(this->elements[ind] * this->elements[ind]);    
04423   }
04424 
04425 
04426   template<class Treal>
04427     void Matrix<Treal>::
04428     frobThreshLowestLevel
04429     (Treal const threshold, Matrix<Treal> * ErrorMatrix) {
04430     if (ErrorMatrix) {
04431       if ((!this->is_zero() && this->frobSquared() <= threshold) ||
04432           (!ErrorMatrix->is_zero() && 
04433            ErrorMatrix->frobSquared() > threshold)) 
04434         Matrix<Treal>::swap(*this,*ErrorMatrix);
04435     }      
04436     else if (!this->is_zero() && this->frobSquared() <= threshold) 
04437       this->clear();
04438   }
04439 
04440   template<class Treal>
04441     void Matrix<Treal>::
04442     frobThreshElementLevel
04443     (Treal const threshold, Matrix<Treal> * ErrorMatrix) {
04444     assert(!this->is_empty());
04445     bool thisMatIsZero = true;
04446     if (ErrorMatrix) {
04447       assert(!ErrorMatrix->is_empty());
04448       bool EMatIsZero = true;
04449       if (!ErrorMatrix->is_zero() || !this->is_zero()) {
04450         if (ErrorMatrix->is_zero())
04451           ErrorMatrix->allocate();
04452         if (this->is_zero())
04453           this->allocate();
04454         for (int ind = 0; ind < this->nElements(); ind++) {
04455           if ( this->elements[ind] != 0 ) {
04456             assert(ErrorMatrix->elements[ind] == 0);
04457             // ok, let's check if we want to move the element to the error matrix
04458             if ( this->elements[ind] * this->elements[ind] <= threshold ) {
04459               // we want to move the element
04460               ErrorMatrix->elements[ind] = this->elements[ind];
04461               this->elements[ind] = 0;
04462               EMatIsZero = false; // at least one element is nonzero
04463             }
04464             else
04465               thisMatIsZero = false; // at least one element is nonzero 
04466             continue;
04467           }
04468           if ( ErrorMatrix->elements[ind] != 0 ) {
04469             assert(this->elements[ind] == 0);
04470             // ok, let's check if we want to move the element from the error matrix
04471             if ( ErrorMatrix->elements[ind] * ErrorMatrix->elements[ind] > threshold ) {
04472               // we want to move the element
04473               this->elements[ind] = ErrorMatrix->elements[ind];
04474               ErrorMatrix->elements[ind] = 0;
04475               thisMatIsZero = false; // at least one element is nonzero 
04476             }
04477             else
04478               EMatIsZero = false; // at least one element is nonzero        
04479           }
04480         }
04481         if (thisMatIsZero) {
04482 #if 0
04483           for (int ind = 0; ind < this->nElements(); ind++) 
04484             assert( this->elements[ind] == 0);
04485 #endif
04486           this->clear();
04487         }
04488         if (EMatIsZero) {
04489 #if 0
04490           for (int ind = 0; ind < this->nElements(); ind++) 
04491             assert( ErrorMatrix->elements[ind] == 0);
04492 #endif
04493           ErrorMatrix->clear();
04494         }
04495       }
04496     }      
04497     else
04498       if (!this->is_zero()) {
04499         for (int ind = 0; ind < this->nElements(); ind++) {
04500           if ( this->elements[ind] * this->elements[ind] <= threshold ) 
04501             /* FIXME BUG? EMANUEL LOOK AT THIS! */
04502             // this->elements[ind] == 0; OLD CODE. Compiler complained that "statement has no effect".
04503             this->elements[ind] = 0;  /* New code. Changed from == to =. */
04504           else
04505             thisMatIsZero = false;
04506         }
04507         if (thisMatIsZero)
04508           this->clear();
04509       }    
04510   }
04511   
04512 
04513   
04514   /********** End of help functions for thresholding */
04515 
04516   /* C = beta * C + alpha * A * B where only the upper triangle of C is */
04517   /* referenced and updated */
04518   template<class Treal>
04519     void Matrix<Treal>:: 
04520     gemm_upper_tr_only(const bool tA, const bool tB, 
04521                        const Treal alpha, 
04522                        const Matrix<Treal>& A, 
04523                        const Matrix<Treal>& B, 
04524                        const Treal beta, 
04525                        Matrix<Treal>& C) {
04526     /* FIXME: Symmetry is not utilized. */
04527     Matrix<Treal>::gemm(tA, tB, alpha, A, B, beta, C);
04528     C.nosymToSym();
04529   }
04530 
04531   /* A = alpha * A * Z or A = alpha * Z * A where A is symmetric, */
04532   /* Z is upper triangular and */
04533   /* only the upper triangle of the result is calculated */
04534   /* side == 'R' for A = alpha * A * Z */
04535   /* side == 'L' for A = alpha * Z * A */
04536   template<class Treal>
04537     void Matrix<Treal>:: 
04538     sytr_upper_tr_only(char const side, const Treal alpha,
04539                        Matrix<Treal>& A,
04540                        const Matrix<Treal>& Z) {
04541     /* FIXME: Symmetry is not utilized. */
04542     A.symToNosym(); 
04543     Matrix<Treal>::trmm(side, 'U', false, alpha, Z, A);
04544     A.nosymToSym();
04545   }
04546   
04547   /* The result B is assumed to be symmetric, i.e. only upper triangle */
04548   /* calculated and hence only upper triangle of input B is used. */
04549   template<class Treal>
04550     void Matrix<Treal>:: 
04551     trmm_upper_tr_only(const char side, const char uplo, 
04552                        const bool tA, const Treal alpha, 
04553                        const Matrix<Treal>& A, 
04554                        Matrix<Treal>& B) {
04555     /* FIXME: Symmetry is not utilized. */
04556     assert(tA);
04557     B.symToNosym();
04558     Matrix<Treal>::trmm(side, uplo, tA, alpha, A, B);
04559     B.nosymToSym();
04560   }
04561 
04562   /* A = Z' * A * Z or A = Z * A * Z' */
04563   /* where A is symmetric and Z is (nonunit) upper triangular */
04564   /* side == 'R' for A = Z' * A * Z */
04565   /* side == 'L' for A = Z * A * Z' */
04566   template<class Treal>
04567     void Matrix<Treal>:: 
04568     trsytriplemm(char const side, 
04569                  const Matrix<Treal>& Z,
04570                  Matrix<Treal>& A) {
04571     if (side == 'R') {
04572       Matrix<Treal>::
04573 	sytr_upper_tr_only('R', 1.0, A, Z);
04574       Matrix<Treal>::
04575 	trmm_upper_tr_only('L', 'U', true, 1.0, Z, A);
04576     }
04577     else {
04578       assert(side == 'L');
04579       Matrix<Treal>::
04580 	sytr_upper_tr_only('L', 1.0, A, Z);
04581       Matrix<Treal>::
04582 	trmm_upper_tr_only('R', 'U', true, 1.0, Z, A);
04583     }
04584   }
04585 
04586 
04587   template<class Treal>
04588     Treal Matrix<Treal>::frob_squared_thresh
04589     (Treal const threshold, Matrix<Treal> * ErrorMatrix) {
04590     assert(!this->is_empty());
04591     if (ErrorMatrix && ErrorMatrix->is_empty()) {
04592       ErrorMatrix->resetRows(this->rows);
04593       ErrorMatrix->resetCols(this->cols);
04594     }
04595     Treal fs = frobSquared();
04596     if (fs < threshold) {
04597       if (ErrorMatrix)
04598         Matrix<Treal>::swap(*this, *ErrorMatrix);
04599       return fs;
04600     }
04601     else
04602       return 0;
04603   }
04604 
04605 
04606   template<class Treal> 
04607     void Matrix<Treal>::
04608     inch(const Matrix<Treal>& A,
04609          Matrix<Treal>& Z,
04610          const Treal threshold, const side looking,
04611          const inchversion version) {
04612     assert(!A.is_empty());
04613     if (A.nrows() != A.ncols()) 
04614       throw Failure("Matrix<Treal>::inch: Matrix must be quadratic!");
04615     if (A.is_zero())
04616       throw Failure("Matrix<Treal>::inch: Matrix is zero! "
04617                     "Not possible to compute inverse cholesky.");
04618     Z = A;
04619     int info;
04620     potrf("U", &A.nrows(), Z.elements, &A.nrows(), &info);
04621     if (info > 0)
04622       throw Failure("Matrix<Treal>::inch: potrf  returned info > 0. The matrix is not positive definite.");
04623     if (info < 0)
04624       throw Failure("Matrix<Treal>::inch: potrf  returned info < 0");
04625       
04626     trtri("U", "N", &A.nrows(), Z.elements, &A.nrows(), &info);
04627     if (info > 0)
04628       throw Failure("Matrix<Treal>::inch: trtri  returned info > 0. The matrix is not positive definite.");
04629     if (info < 0)
04630       throw Failure("Matrix<Treal>::inch: trtri  returned info < 0");
04631     /* Fill lower triangle with zeroes just to be safe */
04632     trifulltofull(Z.elements, A.nrows()); 
04633   }
04634 
04635   template<class Treal>
04636     void Matrix<Treal>:: 
04637     gersgorin(Treal& lmin, Treal& lmax) const {
04638     assert(!this->is_empty());
04639     if (this->nScalarsRows() == this->nScalarsCols()) {
04640       int size = this->nScalarsCols();
04641       Treal* colsums = new Treal[size];
04642       Treal* diag    = new Treal[size];
04643       for (int ind = 0; ind < size; ind++)
04644         colsums[ind] = 0;
04645       this->add_abs_col_sums(colsums);
04646       this->get_diagonal(diag);
04647       Treal tmp1 = colsums[0] - template_blas_fabs(diag[0]);
04648       Treal tmp2;
04649       lmin = diag[0] - tmp1;
04650       lmax = diag[0] + tmp1;
04651       for (int col = 1; col < size; col++) {
04652         tmp1 = colsums[col] - template_blas_fabs(diag[col]);
04653         tmp2 = diag[col] - tmp1;
04654         lmin = (tmp2 < lmin ? tmp2 : lmin);
04655         tmp2 = diag[col] + tmp1;
04656         lmax = (tmp2 > lmax ? tmp2 : lmax);
04657       }
04658       delete[] diag;
04659       delete[] colsums;
04660     }
04661     else
04662       throw Failure("Matrix<Treal>::gersgorin(Treal, Treal): Matrix must be quadratic");
04663   }
04664 
04665 
04666   template<class Treal> 
04667     void Matrix<Treal>::
04668     add_abs_col_sums(Treal* sums) const {
04669     assert(sums);
04670     if (!this->is_zero()) 
04671       for (int col = 0; col < this->ncols(); col++) 
04672         for (int row = 0; row < this->nrows(); row++) 
04673           sums[col] += template_blas_fabs((*this)(row,col));
04674   }
04675  
04676   template<class Treal> 
04677     void Matrix<Treal>::
04678     get_diagonal(Treal* diag) const {
04679     assert(diag);
04680     assert(this->nrows() == this->ncols());
04681     if (this->is_zero()) 
04682       for (int rc = 0; rc < this->nScalarsCols(); rc++)
04683         diag[rc] = 0;
04684     else 
04685       for (int rc = 0; rc < this->ncols(); rc++) 
04686         diag[rc] = (*this)(rc,rc);
04687   }
04688 
04689 
04690     /* New routines above */
04691 
04692 #if 0 /* OLD ROUTINES */   
04693 
04694 
04695   
04696 
04697 
04698 
04699  
04700    
04701 
04702   
04703 
04704 
04705   
04706 
04707 
04708 
04709 
04710 
04711   
04712 
04713   template<class Treal>
04714     void Matrix<Treal>::trim_memory_usage() {
04715     if (this->is_zero() && this->cap > 0) {
04716       delete[] this->elements;
04717       this->elements = NULL;
04718       this->cap = 0;
04719       this->nel = 0;
04720     }
04721     else if (this->cap > this->nel) {
04722       Treal* tmp = new Treal[this->nel];
04723       for (int i = 0; i < this->nel; i++) 
04724         tmp[i] = this->elements[i];
04725       delete[] this->elements;
04726       this->cap = this->nel;
04727       this->elements = tmp;
04728     }
04729     assert(this->cap == this->nel);
04730   }
04731 
04732 
04733 
04734 #if 1
04735 
04736   template<class Treal>
04737     void Matrix<Treal>::
04738     write_to_buffer_count(int& zb_length, int& vb_length) const {
04739     if (this->is_zero()) {
04740       ++zb_length;      
04741     }
04742     else {
04743       ++zb_length;      
04744       vb_length += this->nel;
04745     }
04746   }
04747 
04748   template<class Treal>
04749     void Matrix<Treal>::
04750     write_to_buffer(int* zerobuf, const int zb_length, 
04751                     Treal* valuebuf, const int vb_length,
04752                     int& zb_index, int& vb_index) const {
04753     if (zb_index < zb_length) {
04754       if (this->is_zero()) {
04755         zerobuf[zb_index] = 0;
04756         ++zb_index;     
04757       }
04758       else {
04759         if (vb_index + this->nel < vb_length + 1) {
04760           zerobuf[zb_index] = 1;
04761           ++zb_index;   
04762           for (int i = 0; i < this->nel; i++)
04763             valuebuf[vb_index + i] = this->elements[i];
04764           vb_index += this->nel;
04765         }
04766         else
04767           throw Failure("Matrix<Treal, Telement>::write_to_buffer: "
04768                         "Insufficient space in buffers");   
04769       }
04770     }
04771     else
04772       throw Failure("Matrix<Treal, Telement>::write_to_buffer: "
04773                     "Insufficient space in buffers");   
04774   }
04775 
04776   template<class Treal>
04777     void Matrix<Treal>::
04778     read_from_buffer(int* zerobuf, const int zb_length, 
04779                      Treal* valuebuf, const int vb_length,
04780                      int& zb_index, int& vb_index) {
04781     if (zb_index < zb_length) {
04782       if (zerobuf[zb_index] == 0) {
04783         (*this) = 0;
04784         ++zb_index;     
04785       }
04786       else {
04787         this->content = ful;
04788         this->nel = this->nrows() * this->ncols();
04789         this->assert_alloc();   
04790         if (vb_index + this->nel < vb_length + 1) {
04791           assert(zerobuf[zb_index] == 1);
04792           ++zb_index;   
04793           for (int i = 0; i < this->nel; i++)
04794             this->elements[i] = valuebuf[vb_index + i];
04795           vb_index += this->nel;
04796         }
04797         else
04798           throw Failure("Matrix<Treal, Telement>::read_from_buffer: "
04799                         "Mismatch, buffers does not match matrix");   
04800       }
04801     }
04802     else
04803       throw Failure("Matrix<Treal, Telement>::read_from_buffer: "
04804                     "Mismatch, buffers does not match matrix");   
04805   }
04806   
04807 #endif
04808 
04809 
04810 
04811 
04812 
04813   
04814 
04815 
04816   /* continue here */
04817 
04818 
04819     
04820 
04821     
04822 
04823 
04824 
04825 
04826 
04827 
04828 
04829 
04830 
04831 #if 1
04832 
04833   
04834   
04835 #endif
04836 
04837 #endif     /* END OF OLD ROUTINES */   
04838 
04839 
04840 } /* end namespace mat */
04841 
04842 #endif

Generated on Wed Nov 21 09:32:11 2012 for ergo by  doxygen 1.4.7