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
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 }
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 {
00173 return this->matrixPtr->nnz();
00174 }
00175 inline size_t nvalues() const {
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
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 }
00299 #endif
00300
00301