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