36 #ifndef MAT_MATRIXTRIANGULAR
37 #define MAT_MATRIXTRIANGULAR
56 template<
typename Treal,
typename Tmatrix>
57 class MatrixTriangular :
public MatrixBase<Treal, Tmatrix> {
83 (std::vector<int>
const & rowind,
84 std::vector<int>
const & colind,
85 std::vector<Treal>
const & values,
89 this->
matrixPtr->syAssignFromSparse(rowind, colind, values);
104 (std::vector<int>
const & rowind,
105 std::vector<int>
const & colind,
106 std::vector<Treal>
const & values) {
107 this->
matrixPtr->syAddValues(rowind, colind, values);
112 (std::vector<int>
const & rowind,
113 std::vector<int>
const & colind,
114 std::vector<Treal> & values)
const {
115 this->
matrixPtr->syGetValues(rowind, colind, values);
125 (std::vector<int> & rowind,
126 std::vector<int> & colind,
127 std::vector<Treal> & values)
const {
128 this->
matrixPtr->syGetAllValues(rowind, colind, values);
136 inline void fullmatrix(Treal*
const full,
int const size)
const {
137 this->
matrixPtr->fullmatrix(full, size, size);
142 const Treal threshold,
145 Tmatrix::inch(*SPD.
matrixPtr, *this->matrixPtr,
146 threshold, looking, version);
149 const Treal threshold,
153 Tmatrix::syInch(*SPD.
matrixPtr, *this->matrixPtr,
154 threshold, looking, version);
157 void thresh(Treal
const threshold,
163 Treal
eucl(Treal
const requestedAccuracy,
164 int maxIter = -1)
const;
172 inline size_t nnz()
const {
195 template<
typename TRule>
197 this->
matrixPtr->trSetElementsByRule(rule);
220 template<
typename Treal,
typename Tmatrix>
221 inline MatrixTriangular<Treal, Tmatrix>&
222 MatrixTriangular<Treal, Tmatrix>::operator+=
225 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
230 template<
typename Treal,
typename Tmatrix>
236 this->frob_thresh(threshold);
239 throw Failure(
"MatrixTriangular<Treal, Tmatrix>::"
240 "thresh(Treal const, "
242 "Thresholding not imlpemented for selected norm");
246 template<
typename Treal,
typename Tmatrix>
248 eucl(Treal
const requestedAccuracy,
257 maxIter = this->get_nrows() * 100;
260 lan(ztz, guess, maxIter);
261 lan.setRelTol( 100 * std::numeric_limits<Treal>::epsilon() );
268 if ( euclInt.
low() < 0 )
270 if ( euclInt.
length() / 2.0 > requestedAccuracy ) {
271 std::cout <<
"req: " << requestedAccuracy
272 <<
" obt: " << euclInt.
length() / 2.0 << std::endl;
273 throw std::runtime_error(
"Desired accuracy not obtained in Lanczos.");
280 template<
typename Treal,
typename Tmatrix>
284 return TruncObj.
run( threshold );
289 template<
typename Treal,
typename Tmatrix>
294 return TruncObj.
run( threshold );