1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/libraries/libmbmath/ccmh.h,v 1.40 2017/01/12 14:43:53 masarati Exp $ */ 2 /* 3 * MBDyn (C) is a multibody analysis code. 4 * http://www.mbdyn.org 5 * 6 * Copyright (C) 2003-2017 7 * 8 * This code is a partial merge of HmFe and MBDyn. 9 * 10 * Pierangelo Masarati <masarati@aero.polimi.it> 11 * Paolo Mantegazza <mantegazza@aero.polimi.it> 12 * Marco Morandini <morandini@aero.polimi.it> 13 * 14 * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano 15 * via La Masa, 34 - 20156 Milano, Italy 16 * http://www.aero.polimi.it 17 * 18 * Changing this copyright notice is forbidden. 19 * 20 * This program is free software; you can redistribute it and/or modify 21 * it under the terms of the GNU General Public License as published by 22 * the Free Software Foundation (version 2 of the License). 23 * 24 * 25 * This program is distributed in the hope that it will be useful, 26 * but WITHOUT ANY WARRANTY; without even the implied warranty of 27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 28 * GNU General Public License for more details. 29 * 30 * You should have received a copy of the GNU General Public License 31 * along with this program; if not, write to the Free Software 32 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 33 */ 34 35 #ifndef CColMatrixHandler_hh 36 #define CColMatrixHandler_hh 37 38 #include <vector> 39 40 #include "myassert.h" 41 #include "solman.h" 42 #include "spmh.h" 43 44 /* Sparse Matrix in columns form */ 45 template <int off> 46 class CColMatrixHandler : public CompactSparseMatrixHandler_tpl<off> { 47 private: 48 #ifdef DEBUG IsValid(void)49 void IsValid(void) const { 50 NO_OP; 51 }; 52 #endif /* DEBUG */ 53 54 public: 55 CColMatrixHandler(std::vector<doublereal>& x, 56 const std::vector<integer>& i, 57 const std::vector<integer>& p); 58 59 virtual ~CColMatrixHandler(); 60 61 /* used by MultiThreadDataManager to duplicate the storage array 62 * while preserving the CC indices */ 63 CompactSparseMatrixHandler *Copy(void) const; 64 65 public: operator()66 doublereal & operator()(integer i_row, integer i_col) { 67 ASSERTMSGBREAK(i_row > 0 && i_row <= SparseMatrixHandler::iGetNumRows(), 68 "Error in CColMatrixHandler::operator(), " 69 "row index out of range"); 70 ASSERTMSGBREAK(i_col > 0 && i_col <= SparseMatrixHandler::iGetNumCols(), 71 "Error in CColMatrixHandler::operator(), " 72 "col index out of range"); 73 i_row--; 74 integer row_begin = CompactSparseMatrixHandler_tpl<off>::Ap[i_col - 1] - off; 75 integer row_end = CompactSparseMatrixHandler_tpl<off>::Ap[i_col] - off - 1; 76 integer idx; 77 integer row; 78 79 if (row_begin == (CompactSparseMatrixHandler_tpl<off>::Ap[i_col] - off) 80 || (CompactSparseMatrixHandler_tpl<off>::Ai[row_begin] - off) > i_row 81 || (CompactSparseMatrixHandler_tpl<off>::Ai[row_end] - off) < i_row) 82 { 83 /* matrix must be rebuilt */ 84 throw MatrixHandler::ErrRebuildMatrix(MBDYN_EXCEPT_ARGS); 85 } 86 87 while (row_end >= row_begin) { 88 idx = (row_begin + row_end)/2; 89 row = CompactSparseMatrixHandler_tpl<off>::Ai[idx] - off; 90 if (i_row < row) { 91 row_end = idx - 1; 92 } else if (i_row > row) { 93 row_begin = idx + 1; 94 } else { 95 return CompactSparseMatrixHandler_tpl<off>::Ax[idx]; 96 } 97 } 98 99 /* matrix must be rebuilt */ 100 throw MatrixHandler::ErrRebuildMatrix(MBDYN_EXCEPT_ARGS); 101 }; 102 operator()103 const doublereal& operator () (integer i_row, integer i_col) const { 104 ASSERTMSGBREAK(i_row > 0 && i_row <= SparseMatrixHandler::iGetNumRows(), 105 "Error in CColMatrixHandler::operator(), " 106 "row index out of range"); 107 ASSERTMSGBREAK(i_col > 0 && i_col <= SparseMatrixHandler::iGetNumCols(), 108 "Error in CColMatrixHandler::operator(), " 109 "col index out of range"); 110 i_row--; 111 integer row_begin = CompactSparseMatrixHandler_tpl<off>::Ap[i_col - 1] - off; 112 integer row_end = CompactSparseMatrixHandler_tpl<off>::Ap[i_col] - off - 1; 113 integer idx; 114 integer row; 115 116 if (row_begin == CompactSparseMatrixHandler_tpl<off>::Ap[i_col] - off 117 || CompactSparseMatrixHandler_tpl<off>::Ai[row_begin] - off > i_row 118 || CompactSparseMatrixHandler_tpl<off>::Ai[row_end] - off < i_row) 119 { 120 return ::Zero1; 121 } 122 123 while (row_end >= row_begin) { 124 idx = (row_begin + row_end)/2; 125 row = CompactSparseMatrixHandler_tpl<off>::Ai[idx] - off; 126 if (i_row < row) { 127 row_end = idx - 1; 128 } else if (i_row > row) { 129 row_begin = idx + 1; 130 } else { 131 return CompactSparseMatrixHandler_tpl<off>::Ax[idx]; 132 } 133 } 134 135 return ::Zero1; 136 }; 137 138 void Resize(integer ir, integer ic); 139 140 /* Estrae una colonna da una matrice */ 141 VectorHandler& GetCol(integer icol, VectorHandler& out) const; 142 143 /* Moltiplica per uno scalare e somma a una matrice */ 144 MatrixHandler& MulAndSumWithShift(MatrixHandler& out, 145 doublereal s = 1., 146 integer drow = 0, integer dcol = 0) const; 147 148 MatrixHandler& FakeThirdOrderMulAndSumWithShift(MatrixHandler& out, 149 std::vector<bool> b, doublereal s = 1., 150 integer drow = 0, integer dcol = 0) const; 151 }; 152 153 #endif /* CColMatrixHandler_hh */ 154 155