00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #if !defined(_SPARSE_MATRIX_H_)
00029 #define _SPARSE_MATRIX_H_ 1
00030
00036 #include <stdio.h>
00037 #include <vector>
00038 #include <algorithm>
00039
00040
00041 #include "realtype.h"
00042 #include "matrix_typedefs.h"
00043 #include "basisinfo.h"
00044 #include "sparse_pattern.h"
00045
00046 #if !defined(BEGIN_NAMESPACE)
00047 #define BEGIN_NAMESPACE(x) namespace x {
00048 #define END_NAMESPACE(x) }
00049 #endif
00050
00051 BEGIN_NAMESPACE(Dft)
00052
00053
00054 class SparseMatrix {
00055 class Exception : public std::exception {
00056 const char *msg;
00057 public:
00058 explicit Exception(const char *msg_) : msg(msg_) {}
00059 virtual const char *what() const throw() { return msg; }
00060 };
00061
00062 const SparsePattern& pattern;
00063 ergo_real **columns;
00064 int **offsets;
00065 int **his;
00066 int *cnt;
00067 int n;
00069 void createOffsets(const SparsePattern& pattern);
00070 public:
00073 explicit SparseMatrix(const SparsePattern& pattern_);
00074 SparseMatrix(const SparsePattern& pattern_,
00075 const symmMatrix& m, const int *aoMap,
00076 std::vector<int> const & permutationHML);
00077
00078 ~SparseMatrix() {
00079 for(int i=0; i<n; i++) {
00080 delete [](columns[i]);
00081 delete [](offsets[i]);
00082 delete [](his[i]);
00083 }
00084 delete []columns;
00085 delete []offsets;
00086 delete []his;
00087 delete []cnt;
00088 }
00089
00090 void print(const char *title) const;
00091
00093 void addSymmetrizedTo(symmMatrix& sMat,
00094 const int *aoMap,
00095 std::vector<int> const & permutationHML) const;
00096
00100 void add(int row, int col, ergo_real val) {
00101 ergo_real *columnData = columns[col];
00102 const int *hi = his[col];
00103 int idx;
00104 for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx);
00105
00106 if(idx >= cnt[col])
00107 throw Exception("SparseMatrix::add called with incorrect args");
00108 int offset = offsets[col][idx];
00109
00110
00111
00112 columnData[row-offset] += val;
00113 }
00114
00115
00116
00117
00118
00119 ergo_real at(int row, int col) const {
00120 const ergo_real *columnData = columns[col];
00121 const int *hi = his[col];
00122 int idx; for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx);
00123 if(idx >= cnt[col])
00124 throw Exception("SparseMatrix::at called with incorrect args");
00125
00126 int offset = offsets[col][idx];
00127
00128
00129
00130 return columnData[row-offset];
00131 }
00132 };
00133
00134
00135 END_NAMESPACE(Dft)
00136
00137 void
00138 getrho_blocked_lda(int nbast, const Dft::SparseMatrix& dmat,
00139 const ergo_real * gao,
00140 const int* nblocks, const int (*iblocks)[2],
00141 int ldaib, ergo_real *tmp, int nvclen, ergo_real *rho);
00142 void
00143 getrho_blocked_gga(int nbast, const Dft::SparseMatrix& dmat,
00144 const ergo_real * gao,
00145 const int* nblocks, const int (*iblocks)[2],
00146 int ldaib, ergo_real *tmp, int nvclen,
00147 ergo_real *rho, ergo_real (*grad)[3]);
00148
00149 #endif