1 /*! \file RectBlockMatrix.hpp 2 \brief Block storage for Rectangular matrix 3 \author Atsushi Suzuki, Laboratoire Jacques-Louis Lions 4 \date Jun. 11th 2013 5 \date Jul. 12th 2015 6 \date Nov. 30th 2016 7 */ 8 9 // This file is part of Dissection 10 // 11 // Dissection is free software: you can redistribute it and/or modify 12 // it under the terms of the GNU General Public License as published by 13 // the Free Software Foundation, either version 3 of the License, or 14 // (at your option) any later version. 15 // 16 // Linking Dissection statically or dynamically with other modules is making 17 // a combined work based on Disssection. Thus, the terms and conditions of 18 // the GNU General Public License cover the whole combination. 19 // 20 // As a special exception, the copyright holders of Dissection give you 21 // permission to combine Dissection program with free software programs or 22 // libraries that are released under the GNU LGPL and with independent modules 23 // that communicate with Dissection solely through the Dissection-fortran 24 // interface. You may copy and distribute such a system following the terms of 25 // the GNU GPL for Dissection and the licenses of the other code concerned, 26 // provided that you include the source code of that other code when and as 27 // the GNU GPL requires distribution of source code and provided that you do 28 // not modify the Dissection-fortran interface. 29 // 30 // Note that people who make modified versions of Dissection are not obligated 31 // to grant this special exception for their modified versions; it is their 32 // choice whether to do so. The GNU General Public License gives permission to 33 // release a modified version without this exception; this exception also makes 34 // it possible to release a modified version which carries forward this 35 // exception. If you modify the Dissection-fortran interface, this exception 36 // does not apply to your modified version of Dissection, and you must remove 37 // this exception when you distribute your modified version. 38 // 39 // This exception is an additional permission under section 7 of the GNU 40 // General Public License, version 3 ("GPLv3") 41 // 42 // Dissection is distributed in the hope that it will be useful, 43 // but WITHOUT ANY WARRANTY; without even the implied warranty of 44 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 45 // GNU General Public License for more details. 46 // 47 // You should have received a copy of the GNU General Public License 48 // along with Dissection. If not, see <http://www.gnu.org/licenses/>. 49 // 50 51 #ifndef _ALGEBRA_RECTBLOCKMATRIX_HPP 52 #define _ALGEBRA_RECTBLOCKMATRIX_HPP 53 #include <cstdio> 54 #include <cassert> 55 #include <cmath> 56 #include <vector> 57 58 using std::vector; 59 60 template<typename T> 61 class RectBlockMatrix 62 { 63 public: 64 RectBlockMatrix()65 RectBlockMatrix() 66 { 67 _dim_r = 0; 68 _dim_c = 0; 69 _block_size = 0; 70 _block_size_last_r = 0; 71 _block_size_last_c = 0; 72 _num_blocks_r = 0; 73 _num_blocks_c = 0; 74 _deallocation_status = true; 75 } 76 RectBlockMatrix(int dim_r,int dim_c,int block_size)77 RectBlockMatrix(int dim_r, int dim_c, int block_size) 78 { 79 #ifdef DEBUG_SQUAREBLOCKMATRIX 80 fprintf(stderr, "%s %d : constructor with args %d %d\n", 81 __FILE__, __LINE__, dim, block_size); 82 #endif 83 init(dim_r, dim_c, block_size, 0); 84 allocate(); 85 _deallocation_status = false; 86 } 87 RectBlockMatrix(int dim_r,int dim_c,double * aa)88 RectBlockMatrix(int dim_r, int dim_c, double *aa) 89 { 90 #ifdef DEBUG_SQUAREBLOCKMATRIX 91 fprintf(stderr, "%s %d : constructor with args %d %d\n", 92 __FILE__, __LINE__, dim, block_size); 93 #endif 94 _block_size = dim_r > dim_c ? dim_r : dim_c; 95 _dim_c = dim_c; 96 _dim_r = dim_r; 97 _num_blocks_r = 1; 98 _num_blocks_c = 1; 99 100 _coefs = new T*[1]; 101 _coefs[0] = aa; 102 _block_sizes_r.resize(1); 103 _block_sizes_c.resize(1); 104 _indexblock_r.resize(2); 105 _indexblock_c.resize(2); 106 _block_sizes_r[0] = _dim_r; 107 _block_sizes_c[0] = _dim_c; 108 _indexblock_r[0] = 0; 109 _indexblock_r[1] = _dim_r; 110 _indexblock_c[0] = 0; 111 _indexblock_c[1] = _dim_c; 112 _allocation_status = new bool[1]; 113 _allocation_status[0] = false; // _coefs[0] is not allowed to free 114 _isdecomposed_c = false; 115 _deallocation_status = true; 116 } 117 118 void init(int dim_r, int dim_c, int block_size, int first_block = 0); 119 void allocateBlock(int i, int j); 120 void allocate(); 121 void free(int i, int j); 122 void free(); 123 ~RectBlockMatrix()124 ~RectBlockMatrix() { 125 free(); 126 if (!_deallocation_status) { 127 delete [] _coefs; 128 delete [] _allocation_status; 129 } 130 _deallocation_status = true; 131 _block_sizes_r.clear(); 132 _block_sizes_c.clear(); 133 _num_blocks_r = 0; // ? to avoid double free 134 _num_blocks_c = 0; // ? to avoid double free 135 _dim_r = 0; 136 _dim_c = 0; 137 } 138 dimension_r() const139 int dimension_r() const { return (int)_dim_r; } dimension_c() const140 int dimension_c() const { return (int)_dim_c; } block_size() const141 int block_size() const { return (int)_block_size; } num_blocks_r() const142 int num_blocks_r() const { return _num_blocks_r; } num_blocks_c() const143 int num_blocks_c() const { return _num_blocks_c; } nbRows() const144 int nbRows() const { return _dim_r; } nbColumns() const145 int nbColumns() const { return _dim_c; } 146 block_size_last_r() const147 int block_size_last_r() const { return (int)_block_size_last_r; } block_size_last_c() const148 int block_size_last_c() const { return (int)_block_size_last_c; } 149 nrowBlock(int i) const150 int nrowBlock(int i) const { return _block_sizes_r[i]; } ncolBlock(int i) const151 int ncolBlock(int i) const { return _block_sizes_c[i]; } 152 IndexBlock_r(int i) const153 int IndexBlock_r(int i) const { return _indexblock_r[i]; } IndexBlock_c(int i) const154 int IndexBlock_c(int i) const { return _indexblock_c[i]; } 155 T* addrCoefBlock(int i, int j); 156 157 void ZeroClear(); 158 159 T& operator () (int i, int j); 160 161 const T& operator () (int i, int j) const; 162 clone() const163 RectBlockMatrix<T>* clone() const 164 { 165 RectBlockMatrix<T> *ret = new RectBlockMatrix<T>; 166 ret->copy(*this); 167 return(ret); 168 } 169 copy(const RectBlockMatrix & A)170 void copy ( const RectBlockMatrix& A ) 171 { 172 _block_size = A._block_size; 173 _block_size2 = A._block_size2; 174 _num_blocks = A._num_blocks; 175 _block_sizes_r = A._block_sizes_r; 176 _block_sizes_c = A._block_sizes_c; 177 _dim_r = A._dim_r; 178 _dim_c = A._dim_c; 179 _block_size_last_r = A._block_size_last_r; 180 _block_size_last_c = A._block_size_last_c; 181 _num_blocks_r = A._num_blocks_r; 182 _num_blocks_c = A._num_blocks_c; 183 _coefs = A._coefs; 184 _allocation_status = A._allocation_status; 185 _deallocation_status = A._deallocation_status; 186 } 187 188 private: 189 int _block_size; 190 int _block_size2; 191 int _num_blocks; 192 vector<int> _block_sizes_r; 193 vector<int> _block_sizes_c; 194 vector<int> _indexblock_r; 195 vector<int> _indexblock_c; 196 int _dim_r; 197 int _dim_c; 198 int _block_size_last_r; 199 int _block_size_last_c; 200 int _block_size_last0_c; 201 int _block_size_last1_c; 202 int _num_blocks_r; 203 int _num_blocks_c; 204 int _num_blocks0_c; 205 int _num_blocks1_c; 206 int _dim0_c; 207 int _dim1_c; 208 T** _coefs; 209 bool *_allocation_status; 210 bool _deallocation_status; 211 bool _isdecomposed_c; 212 }; // End class RectMatrix 213 214 #endif 215