1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
|
/* Ergo, version 3.8.2, a program for linear scaling electronic structure
* calculations.
* Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
* and Anastasia Kruchinina.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Primary academic reference:
* Ergo: An open-source program for linear-scaling electronic structure
* calculations,
* Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
* Kruchinina,
* SoftwareX 7, 107 (2018),
* <http://dx.doi.org/10.1016/j.softx.2018.03.005>
*
* For further information about Ergo, see <http://www.ergoscf.org>.
*/
/** @file MatrixTriangular.h Triangular matrix class
*
* Copyright(c) Emanuel Rubensson 2006
*
* @author Emanuel Rubensson @a responsible @a author
* @date January 2006
*
*/
#ifndef MAT_MATRIXTRIANGULAR
#define MAT_MATRIXTRIANGULAR
#include <stdexcept>
#include "MatrixBase.h"
namespace mat {
/** Upper non-unit triangular matrix
*
*
* This class belongs to the matrix API
*
* Treal: Type for real numbers
*
* Tmatrix: The matrix class
*
* @see MatrixBase
* @see MatrixGeneral
* @see MatrixSymmetric
*
*
*/
template<typename Treal, typename Tmatrix>
class MatrixTriangular : public MatrixBase<Treal, Tmatrix> {
public:
typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType;
typedef Treal real;
MatrixTriangular()
:MatrixBase<Treal, Tmatrix>() {} /**< Default constructor */
explicit MatrixTriangular(const MatrixTriangular<Treal, Tmatrix>& tri)
:MatrixBase<Treal, Tmatrix>(tri) {} /**< Copy constructor */
MatrixTriangular<Treal, Tmatrix>&
operator=(const MatrixTriangular<Treal, Tmatrix>& tri) {
MatrixBase<Treal, Tmatrix>::operator=(tri);
return *this;
}
inline MatrixTriangular<Treal, Tmatrix>& operator=(int const k) {
*this->matrixPtr = k;
return *this;
} /**< Set matrix to zero or identity: A = 0 or A = 1
*
* Only zero and one are valid arguments.
*
*
*/
inline void assign_from_sparse
(std::vector<int> const & rowind,
std::vector<int> const & colind,
std::vector<Treal> const & values,
SizesAndBlocks const & newRows,
SizesAndBlocks const & newCols) {
this->resetSizesAndBlocks(newRows, newCols);
this->matrixPtr->syAssignFromSparse(rowind, colind, values);
}
/**< Assign from sparse matrix given by three vectors.
* The vectors contain row indices, column indices and values.
* The indices start at zero.
* The elements to be added must be in upper triangle.
* Information about sizes and blocks for rows as well as columns
* must also be given.
* @warning All indexing start at zero.
*/
/** Add given set of values to the matrix (+=).
* The values should be in upper triangle.
*/
inline void add_values
(std::vector<int> const & rowind,
std::vector<int> const & colind,
std::vector<Treal> const & values) {
this->matrixPtr->syAddValues(rowind, colind, values);
}
inline void get_values
(std::vector<int> const & rowind,
std::vector<int> const & colind,
std::vector<Treal> & values) const {
this->matrixPtr->syGetValues(rowind, colind, values);
}
/**< Get values given by row and column index lists.
* Input arrays contain row and column indices.
* The wanted elements must in upper triangle.
* The output array contains values for the given indices.
* @warning All indexing start at zero.
*/
inline void get_all_values
(std::vector<int> & rowind,
std::vector<int> & colind,
std::vector<Treal> & values) const {
rowind.resize(0);
colind.resize(0);
values.resize(0);
rowind.reserve(nnz());
colind.reserve(nnz());
values.reserve(nnz());
this->matrixPtr->syGetAllValues(rowind, colind, values);
}
/**< Get all values and corresponding row and column index lists,
* in matrix. Only upper triangle values are returned.
* @warning All indexing start at zero.
*/
#if 0
inline void fullmatrix(Treal* const full, int const size) const {
this->matrixPtr->fullmatrix(full, size, size);
} /* FIXME? Should triangular matrix always have zeros below diagonal? */
#endif
inline void inch(const MatrixGeneral<Treal, Tmatrix>& SPD,
const Treal threshold,
const side looking = left,
const inchversion version = unstable) {
Tmatrix::inch(*SPD.matrixPtr, *this->matrixPtr,
threshold, looking, version);
}
inline void inch(const MatrixSymmetric<Treal, Tmatrix>& SPD,
const Treal threshold,
const side looking = left,
const inchversion version = unstable) {
this->matrixPtr.haveDataStructureSet(true);
Tmatrix::syInch(*SPD.matrixPtr, *this->matrixPtr,
threshold, looking, version);
}
void thresh(Treal const threshold,
normType const norm);
inline Treal frob() const {
return this->matrixPtr->frob();
}
Treal eucl(Treal const requestedAccuracy,
int maxIter = -1) const;
Treal eucl_thresh(Treal const threshold);
Treal eucl_thresh_congr_trans_measure(Treal const threshold,
MatrixSymmetric<Treal, Tmatrix> & trA);
inline void frob_thresh(Treal threshold) {
this->matrixPtr->frob_thresh(threshold);
}
inline size_t nnz() const { /* Note: size_t instead of int here to avoid integer overflow. */
return this->matrixPtr->nnz();
}
inline size_t nvalues() const { /* Note: size_t instead of int here to avoid integer overflow. */
return this->matrixPtr->nvalues();
}
inline void write_to_buffer(void* buffer, const int n_bytes) const {
this->write_to_buffer_base(buffer, n_bytes, matrix_triang);
}
inline void read_from_buffer(void* buffer, const int n_bytes) {
this->read_from_buffer_base(buffer, n_bytes, matrix_triang);
}
void random() {
this->matrixPtr->syRandom();
}
/** Uses rule depending on the row and column indexes to set matrix elements
* The Trule class should have the function "Treal = set(int row,int col)"
* which is used to set the elements.
*/
template<typename TRule>
void setElementsByRule(TRule & rule) {
this->matrixPtr->trSetElementsByRule(rule);
return;
}
/** B += alpha * A */
MatrixTriangular<Treal, Tmatrix>& operator+=
(XY<Treal, MatrixTriangular<Treal, Tmatrix> > const & sm);
std::string obj_type_id() const {return "MatrixTriangular";}
protected:
inline void writeToFileProt(std::ofstream & file) const {
this->writeToFileBase(file, matrix_triang);
}
inline void readFromFileProt(std::ifstream & file) {
this->readFromFileBase(file, matrix_triang);
}
private:
};
/* B += alpha * A */
template<typename Treal, typename Tmatrix>
inline MatrixTriangular<Treal, Tmatrix>&
MatrixTriangular<Treal, Tmatrix>::operator+=
(XY<Treal, MatrixTriangular<Treal, Tmatrix> > const & sm) {
assert(!sm.tB);
Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
return *this;
}
template<typename Treal, typename Tmatrix>
void MatrixTriangular<Treal, Tmatrix>::
thresh(Treal const threshold,
normType const norm) {
switch (norm) {
case frobNorm:
this->frob_thresh(threshold);
break;
default:
throw Failure("MatrixTriangular<Treal, Tmatrix>::"
"thresh(Treal const, "
"normType const): "
"Thresholding not imlpemented for selected norm");
}
}
template<typename Treal, typename Tmatrix>
Treal MatrixTriangular<Treal, Tmatrix>::
eucl(Treal const requestedAccuracy,
int maxIter) const {
VectorType guess;
SizesAndBlocks cols;
this->getCols(cols);
guess.resetSizesAndBlocks(cols);
guess.rand();
mat::ATAMatrix<MatrixTriangular<Treal, Tmatrix>, Treal> ztz(*this);
if (maxIter < 0)
maxIter = this->get_nrows() * 100;
arn::LanczosLargestMagnitudeEig
<Treal, ATAMatrix<MatrixTriangular<Treal, Tmatrix>, Treal>, VectorType>
lan(ztz, guess, maxIter);
lan.setRelTol( 100 * mat::getMachineEpsilon<Treal>() );
lan.run();
Treal eVal = 0;
Treal acc = 0;
lan.getLargestMagnitudeEig(eVal, acc);
Interval<Treal> euclInt( template_blas_sqrt(eVal-acc),
template_blas_sqrt(eVal+acc) );
if ( euclInt.low() < 0 )
euclInt = Interval<Treal>( 0, template_blas_sqrt(eVal+acc) );
if ( euclInt.length() / 2.0 > requestedAccuracy ) {
std::cout << "req: " << (double)requestedAccuracy
<< " obt: " << (double)(euclInt.length() / 2.0) << std::endl;
throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
}
return euclInt.midPoint();
}
#if 1
template<typename Treal, typename Tmatrix>
Treal MatrixTriangular<Treal, Tmatrix>::
eucl_thresh(Treal const threshold) {
EuclTruncationGeneral<MatrixTriangular<Treal, Tmatrix>, Treal> TruncObj( *this );
return TruncObj.run( threshold );
}
#endif
template<typename Treal, typename Tmatrix>
Treal MatrixTriangular<Treal, Tmatrix>::
eucl_thresh_congr_trans_measure(Treal const threshold,
MatrixSymmetric<Treal, Tmatrix> & trA) {
EuclTruncationCongrTransMeasure<MatrixTriangular<Treal, Tmatrix>, MatrixSymmetric<Treal, Tmatrix>, Treal> TruncObj(*this, trA);
return TruncObj.run( threshold );
}
} /* end namespace mat */
#endif
|