MatrixTriangular.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 
00036 #ifndef MAT_MATRIXTRIANGULAR
00037 #define MAT_MATRIXTRIANGULAR
00038 #include <stdexcept>
00039 #include "MatrixBase.h"
00040 namespace mat {
00056   template<typename Treal, typename Tmatrix>
00057     class MatrixTriangular : public MatrixBase<Treal, Tmatrix> {
00058   public:
00059     typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType;
00060 
00061     MatrixTriangular()
00062       :MatrixBase<Treal, Tmatrix>() {} 
00063     explicit MatrixTriangular(const MatrixTriangular<Treal, Tmatrix>& tri)
00064       :MatrixBase<Treal, Tmatrix>(tri) {} 
00066     MatrixTriangular<Treal, Tmatrix>& 
00067       operator=(const MatrixTriangular<Treal, Tmatrix>& tri) {
00068       MatrixBase<Treal, Tmatrix>::operator=(tri);
00069       return *this;
00070     } 
00071     
00072     inline MatrixTriangular<Treal, Tmatrix>& operator=(int const k) {
00073       *this->matrixPtr = k;
00074       return *this;
00075     } 
00082     inline void assign_from_sparse
00083       (std::vector<int> const & rowind, 
00084        std::vector<int> const & colind, 
00085        std::vector<Treal> const & values, 
00086        SizesAndBlocks const & newRows,
00087        SizesAndBlocks const & newCols) {
00088       this->resetSizesAndBlocks(newRows, newCols);
00089       this->matrixPtr->syAssignFromSparse(rowind, colind, values);
00090     }
00103     inline void add_values
00104       (std::vector<int> const & rowind, 
00105        std::vector<int> const & colind, 
00106        std::vector<Treal> const & values) {
00107       this->matrixPtr->syAddValues(rowind, colind, values);
00108     }
00109 
00110 
00111     inline void get_values
00112       (std::vector<int> const & rowind, 
00113        std::vector<int> const & colind, 
00114        std::vector<Treal> & values) const {
00115       this->matrixPtr->syGetValues(rowind, colind, values);
00116     }
00124     inline void get_all_values
00125       (std::vector<int> & rowind, 
00126        std::vector<int> & colind, 
00127        std::vector<Treal> & values) const {
00128       this->matrixPtr->syGetAllValues(rowind, colind, values);
00129     }
00135 #if 0
00136     inline void fullmatrix(Treal* const full, int const size) const {
00137       this->matrixPtr->fullmatrix(full, size, size);
00138     } /* FIXME? Should triangular matrix always have zeros below diagonal? */
00139 #endif
00140 
00141     inline void inch(const MatrixGeneral<Treal, Tmatrix>& SPD, 
00142                      const Treal threshold,
00143                      const side looking = left,
00144                      const inchversion version = unstable) {
00145       Tmatrix::inch(*SPD.matrixPtr, *this->matrixPtr, 
00146                     threshold, looking, version);
00147     }
00148     inline void inch(const MatrixSymmetric<Treal, Tmatrix>& SPD, 
00149                      const Treal threshold,
00150                      const side looking = left,
00151                      const inchversion version = unstable) {
00152       this->matrixPtr.haveDataStructureSet(true);
00153       Tmatrix::syInch(*SPD.matrixPtr, *this->matrixPtr, 
00154                       threshold, looking, version);
00155     }
00156     
00157     void thresh(Treal const threshold,
00158                 normType const norm);
00159 
00160     inline Treal frob() const {
00161       return this->matrixPtr->frob();
00162     }
00163     Treal eucl(Treal const requestedAccuracy,
00164                int maxIter = -1) const;
00165 
00166     Treal eucl_thresh(Treal const threshold);
00167     Treal eucl_thresh_congr_trans_measure(Treal const threshold,
00168                                           MatrixSymmetric<Treal, Tmatrix> & trA);
00169     inline void frob_thresh(Treal threshold) {
00170       this->matrixPtr->frob_thresh(threshold);
00171     }    
00172     inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
00173       return this->matrixPtr->nnz();
00174     }
00175     inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
00176       return this->matrixPtr->nvalues();
00177     }
00178 
00179 
00180     inline void write_to_buffer(void* buffer, const int n_bytes) const {
00181       this->write_to_buffer_base(buffer, n_bytes, matrix_triang);
00182     }
00183     inline void read_from_buffer(void* buffer, const int n_bytes) {
00184       this->read_from_buffer_base(buffer, n_bytes, matrix_triang);
00185     }
00186 
00187     void random() {
00188       this->matrixPtr->syRandom();
00189     }
00190 
00195     template<typename TRule>
00196       void setElementsByRule(TRule & rule) {
00197       this->matrixPtr->trSetElementsByRule(rule);
00198       return;
00199     }
00200 
00202     MatrixTriangular<Treal, Tmatrix>& operator+= 
00203       (XY<Treal, MatrixTriangular<Treal, Tmatrix> > const & sm);
00204  
00205 
00206     std::string obj_type_id() const {return "MatrixTriangular";}
00207     protected:
00208     inline void writeToFileProt(std::ofstream & file) const {
00209       this->writeToFileBase(file, matrix_triang);
00210     }
00211     inline void readFromFileProt(std::ifstream & file) {
00212       this->readFromFileBase(file, matrix_triang);
00213     }
00214 
00215     private:
00216 
00217   };
00218 
00219   /* B += alpha * A */
00220   template<typename Treal, typename Tmatrix>
00221     inline MatrixTriangular<Treal, Tmatrix>& 
00222     MatrixTriangular<Treal, Tmatrix>::operator+= 
00223     (XY<Treal, MatrixTriangular<Treal, Tmatrix> > const & sm) {
00224     assert(!sm.tB);
00225     Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
00226     return *this;
00227   }
00228 
00229 
00230   template<typename Treal, typename Tmatrix>
00231     void MatrixTriangular<Treal, Tmatrix>::
00232     thresh(Treal const threshold,
00233            normType const norm) {
00234     switch (norm) {
00235     case frobNorm:
00236       this->frob_thresh(threshold);
00237       break;
00238     default:
00239       throw Failure("MatrixTriangular<Treal, Tmatrix>::"
00240                     "thresh(Treal const, "
00241                     "normType const): "
00242                     "Thresholding not imlpemented for selected norm");
00243     }
00244   }
00245 
00246   template<typename Treal, typename Tmatrix>
00247     Treal MatrixTriangular<Treal, Tmatrix>::
00248     eucl(Treal const requestedAccuracy,
00249          int maxIter) const {
00250     VectorType guess;
00251     SizesAndBlocks cols;
00252     this->getCols(cols);
00253     guess.resetSizesAndBlocks(cols);
00254     guess.rand();
00255     mat::ATAMatrix<MatrixTriangular<Treal, Tmatrix>, Treal> ztz(*this);    
00256     if (maxIter < 0)
00257       maxIter = this->get_nrows() * 100;
00258     arn::LanczosLargestMagnitudeEig
00259       <Treal, ATAMatrix<MatrixTriangular<Treal, Tmatrix>, Treal>, VectorType>
00260       lan(ztz, guess, maxIter);
00261     lan.setRelTol( 100 * std::numeric_limits<Treal>::epsilon() );
00262     lan.run();
00263     Treal eVal = 0;
00264     Treal acc = 0;
00265     lan.getLargestMagnitudeEig(eVal, acc);
00266     Interval<Treal> euclInt( template_blas_sqrt(eVal-acc),
00267                              template_blas_sqrt(eVal+acc) );    
00268     if ( euclInt.low() < 0 )
00269       euclInt = Interval<Treal>( 0, template_blas_sqrt(eVal+acc) );
00270     if ( euclInt.length() / 2.0 > requestedAccuracy ) {
00271       std::cout << "req: " << requestedAccuracy
00272                 << "  obt: " << euclInt.length() / 2.0 << std::endl;
00273       throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
00274     }
00275     return euclInt.midPoint();
00276   }
00277 
00278 #if 1
00279 
00280   template<typename Treal, typename Tmatrix>
00281     Treal MatrixTriangular<Treal, Tmatrix>::
00282     eucl_thresh(Treal const threshold) {
00283     EuclTruncationGeneral<MatrixTriangular<Treal, Tmatrix>, Treal> TruncObj( *this );
00284     return TruncObj.run( threshold );    
00285   }
00286   
00287 #endif  
00288 
00289   template<typename Treal, typename Tmatrix>
00290     Treal MatrixTriangular<Treal, Tmatrix>::
00291     eucl_thresh_congr_trans_measure(Treal const threshold,
00292                                MatrixSymmetric<Treal, Tmatrix> & trA) {
00293     EuclTruncationCongrTransMeasure<MatrixTriangular<Treal, Tmatrix>, MatrixSymmetric<Treal, Tmatrix>, Treal> TruncObj(*this, trA);
00294     return TruncObj.run( threshold );    
00295   }
00296   
00297 
00298 } /* end namespace mat */
00299 #endif
00300 
00301 

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