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
00038 #ifndef MAT_MATRIXTRIANGULAR
00039 #define MAT_MATRIXTRIANGULAR
00040 #include <stdexcept>
00041 #include "MatrixBase.h"
00042 namespace mat {
00058 template<typename Treal, typename Tmatrix>
00059 class MatrixTriangular : public MatrixBase<Treal, Tmatrix> {
00060 public:
00061 typedef VectorGeneral<Treal, typename Tmatrix::VectorType> VectorType;
00062 typedef Treal real;
00063
00064 MatrixTriangular()
00065 :MatrixBase<Treal, Tmatrix>() {}
00066 explicit MatrixTriangular(const MatrixTriangular<Treal, Tmatrix>& tri)
00067 :MatrixBase<Treal, Tmatrix>(tri) {}
00069 MatrixTriangular<Treal, Tmatrix>&
00070 operator=(const MatrixTriangular<Treal, Tmatrix>& tri) {
00071 MatrixBase<Treal, Tmatrix>::operator=(tri);
00072 return *this;
00073 }
00074
00075 inline MatrixTriangular<Treal, Tmatrix>& operator=(int const k) {
00076 *this->matrixPtr = k;
00077 return *this;
00078 }
00085 inline void assign_from_sparse
00086 (std::vector<int> const & rowind,
00087 std::vector<int> const & colind,
00088 std::vector<Treal> const & values,
00089 SizesAndBlocks const & newRows,
00090 SizesAndBlocks const & newCols) {
00091 this->resetSizesAndBlocks(newRows, newCols);
00092 this->matrixPtr->syAssignFromSparse(rowind, colind, values);
00093 }
00106 inline void add_values
00107 (std::vector<int> const & rowind,
00108 std::vector<int> const & colind,
00109 std::vector<Treal> const & values) {
00110 this->matrixPtr->syAddValues(rowind, colind, values);
00111 }
00112
00113
00114 inline void get_values
00115 (std::vector<int> const & rowind,
00116 std::vector<int> const & colind,
00117 std::vector<Treal> & values) const {
00118 this->matrixPtr->syGetValues(rowind, colind, values);
00119 }
00127 inline void get_all_values
00128 (std::vector<int> & rowind,
00129 std::vector<int> & colind,
00130 std::vector<Treal> & values) const {
00131 rowind.resize(0);
00132 colind.resize(0);
00133 values.resize(0);
00134 rowind.reserve(nnz());
00135 colind.reserve(nnz());
00136 values.reserve(nnz());
00137 this->matrixPtr->syGetAllValues(rowind, colind, values);
00138 }
00144 #if 0
00145 inline void fullmatrix(Treal* const full, int const size) const {
00146 this->matrixPtr->fullmatrix(full, size, size);
00147 }
00148 #endif
00149
00150 inline void inch(const MatrixGeneral<Treal, Tmatrix>& SPD,
00151 const Treal threshold,
00152 const side looking = left,
00153 const inchversion version = unstable) {
00154 Tmatrix::inch(*SPD.matrixPtr, *this->matrixPtr,
00155 threshold, looking, version);
00156 }
00157 inline void inch(const MatrixSymmetric<Treal, Tmatrix>& SPD,
00158 const Treal threshold,
00159 const side looking = left,
00160 const inchversion version = unstable) {
00161 this->matrixPtr.haveDataStructureSet(true);
00162 Tmatrix::syInch(*SPD.matrixPtr, *this->matrixPtr,
00163 threshold, looking, version);
00164 }
00165
00166 void thresh(Treal const threshold,
00167 normType const norm);
00168
00169 inline Treal frob() const {
00170 return this->matrixPtr->frob();
00171 }
00172 Treal eucl(Treal const requestedAccuracy,
00173 int maxIter = -1) const;
00174
00175 Treal eucl_thresh(Treal const threshold);
00176 Treal eucl_thresh_congr_trans_measure(Treal const threshold,
00177 MatrixSymmetric<Treal, Tmatrix> & trA);
00178 inline void frob_thresh(Treal threshold) {
00179 this->matrixPtr->frob_thresh(threshold);
00180 }
00181 inline size_t nnz() const {
00182 return this->matrixPtr->nnz();
00183 }
00184 inline size_t nvalues() const {
00185 return this->matrixPtr->nvalues();
00186 }
00187
00188
00189 inline void write_to_buffer(void* buffer, const int n_bytes) const {
00190 this->write_to_buffer_base(buffer, n_bytes, matrix_triang);
00191 }
00192 inline void read_from_buffer(void* buffer, const int n_bytes) {
00193 this->read_from_buffer_base(buffer, n_bytes, matrix_triang);
00194 }
00195
00196 void random() {
00197 this->matrixPtr->syRandom();
00198 }
00199
00204 template<typename TRule>
00205 void setElementsByRule(TRule & rule) {
00206 this->matrixPtr->trSetElementsByRule(rule);
00207 return;
00208 }
00209
00211 MatrixTriangular<Treal, Tmatrix>& operator+=
00212 (XY<Treal, MatrixTriangular<Treal, Tmatrix> > const & sm);
00213
00214
00215 std::string obj_type_id() const {return "MatrixTriangular";}
00216 protected:
00217 inline void writeToFileProt(std::ofstream & file) const {
00218 this->writeToFileBase(file, matrix_triang);
00219 }
00220 inline void readFromFileProt(std::ifstream & file) {
00221 this->readFromFileBase(file, matrix_triang);
00222 }
00223
00224 private:
00225
00226 };
00227
00228
00229 template<typename Treal, typename Tmatrix>
00230 inline MatrixTriangular<Treal, Tmatrix>&
00231 MatrixTriangular<Treal, Tmatrix>::operator+=
00232 (XY<Treal, MatrixTriangular<Treal, Tmatrix> > const & sm) {
00233 assert(!sm.tB);
00234 Tmatrix::add(sm.A, *sm.B.matrixPtr, *this->matrixPtr);
00235 return *this;
00236 }
00237
00238
00239 template<typename Treal, typename Tmatrix>
00240 void MatrixTriangular<Treal, Tmatrix>::
00241 thresh(Treal const threshold,
00242 normType const norm) {
00243 switch (norm) {
00244 case frobNorm:
00245 this->frob_thresh(threshold);
00246 break;
00247 default:
00248 throw Failure("MatrixTriangular<Treal, Tmatrix>::"
00249 "thresh(Treal const, "
00250 "normType const): "
00251 "Thresholding not imlpemented for selected norm");
00252 }
00253 }
00254
00255 template<typename Treal, typename Tmatrix>
00256 Treal MatrixTriangular<Treal, Tmatrix>::
00257 eucl(Treal const requestedAccuracy,
00258 int maxIter) const {
00259 VectorType guess;
00260 SizesAndBlocks cols;
00261 this->getCols(cols);
00262 guess.resetSizesAndBlocks(cols);
00263 guess.rand();
00264 mat::ATAMatrix<MatrixTriangular<Treal, Tmatrix>, Treal> ztz(*this);
00265 if (maxIter < 0)
00266 maxIter = this->get_nrows() * 100;
00267 arn::LanczosLargestMagnitudeEig
00268 <Treal, ATAMatrix<MatrixTriangular<Treal, Tmatrix>, Treal>, VectorType>
00269 lan(ztz, guess, maxIter);
00270 lan.setRelTol( 100 * mat::getMachineEpsilon<Treal>() );
00271 lan.run();
00272 Treal eVal = 0;
00273 Treal acc = 0;
00274 lan.getLargestMagnitudeEig(eVal, acc);
00275 Interval<Treal> euclInt( template_blas_sqrt(eVal-acc),
00276 template_blas_sqrt(eVal+acc) );
00277 if ( euclInt.low() < 0 )
00278 euclInt = Interval<Treal>( 0, template_blas_sqrt(eVal+acc) );
00279 if ( euclInt.length() / 2.0 > requestedAccuracy ) {
00280 std::cout << "req: " << (double)requestedAccuracy
00281 << " obt: " << (double)(euclInt.length() / 2.0) << std::endl;
00282 throw std::runtime_error("Desired accuracy not obtained in Lanczos.");
00283 }
00284 return euclInt.midPoint();
00285 }
00286
00287 #if 1
00288
00289 template<typename Treal, typename Tmatrix>
00290 Treal MatrixTriangular<Treal, Tmatrix>::
00291 eucl_thresh(Treal const threshold) {
00292 EuclTruncationGeneral<MatrixTriangular<Treal, Tmatrix>, Treal> TruncObj( *this );
00293 return TruncObj.run( threshold );
00294 }
00295
00296 #endif
00297
00298 template<typename Treal, typename Tmatrix>
00299 Treal MatrixTriangular<Treal, Tmatrix>::
00300 eucl_thresh_congr_trans_measure(Treal const threshold,
00301 MatrixSymmetric<Treal, Tmatrix> & trA) {
00302 EuclTruncationCongrTransMeasure<MatrixTriangular<Treal, Tmatrix>, MatrixSymmetric<Treal, Tmatrix>, Treal> TruncObj(*this, trA);
00303 return TruncObj.run( threshold );
00304 }
00305
00306
00307 }
00308 #endif
00309
00310