00001
00003 inline zgbmatrix& zgbmatrix::operator=(const _zgbmatrix& mat)
00004 {
00005 #ifdef CPPL_VERBOSE
00006 std::cerr << "# [MARK] zgbmatrix::operator=(const _zgbmatrix&)"
00007 << std::endl;
00008 #endif//CPPL_VERBOSE
00009
00010 shallow_copy(mat);
00011 return *this;
00012 }
00013
00017
00018
00021 inline zgbmatrix& zgbmatrix::operator+=(const _zgbmatrix& mat)
00022 {
00023 #ifdef CPPL_VERBOSE
00024 std::cerr << "# [MARK] zgbmatrix::operator+=(const _zgbmatrix&)"
00025 << std::endl;
00026 #endif//CPPL_VERBOSE
00027
00028 #ifdef CPPL_DEBUG
00029 if(N!=mat.N || M!=mat.M){
00030 std::cerr << "[ERROR] zgbmatrix::operator+=(_zgbmatrix&)" << std::endl
00031 << "These two matrises can not make a summation." << std::endl
00032 << "Your input was"
00033 << "(" << M <<"x"<< N <<","<< KL <<":"<< KU << ") "<< "+="
00034 << "("<< mat.M <<"x"<< mat.N <<","<< mat.KL <<":"<< mat.KU <<") "
00035 << std::endl;
00036 exit(1);
00037 }
00038 #endif//CPPL_DEBUG
00039
00040 if(KL>=mat.KL && KU>=mat.KU){
00041 for(long i=0; i<M; i++){
00042 for(long j=max(0,i-mat.KL); j<min(N,i+mat.KU+1); j++){
00043 operator()(i,j)+=mat(i,j);
00044 }
00045 }
00046
00047 mat.destroy();
00048 return *this;
00049 }
00050 else{
00051 zgbmatrix newmat(M,N,max(KL,mat.KL),max(KU,mat.KU));
00052 newmat.zero();
00053 for(long i=0; i<M; i++){
00054 for(long j=max(0,i-KL); j<min(N,i+KU+1); j++){
00055 newmat(i,j)+=operator()(i,j);
00056 }
00057 for(long j=max(0,i-mat.KL); j<min(mat.N,i+mat.KU+1); j++){
00058 newmat(i,j)+=mat(i,j);
00059 }
00060 }
00061
00062 swap(*this,newmat);
00063 mat.destroy();
00064 return *this;
00065 }
00066 }
00067
00068
00071 inline zgbmatrix& zgbmatrix::operator-=(const _zgbmatrix& mat)
00072 {
00073 #ifdef CPPL_VERBOSE
00074 std::cerr << "# [MARK] zgbmatrix::operator-=(const _zgbmatrix&)"
00075 << std::endl;
00076 #endif//CPPL_VERBOSE
00077
00078 #ifdef CPPL_DEBUG
00079 if(N!=mat.N || M!=mat.M){
00080 std::cerr << "[ERROR] zgbmatrix::operator-=(_zgbmatrix&)" << std::endl
00081 << "These two matrises can not make a subtraction." << std::endl
00082 << "Your input was"
00083 << "(" << M <<"x"<< N <<","<< KL <<":"<< KU << ") "<< "-="
00084 << "("<< mat.M <<"x"<< mat.N <<","<< mat.KL <<":"<< mat.KU <<") "
00085 << std::endl;
00086 exit(1);
00087 }
00088 #endif//CPPL_DEBUG
00089
00090 if(KL>=mat.KL && KU>=mat.KU){
00091 for(long i=0; i<M; i++){
00092 for(long j=max(0,i-mat.KL); j<min(N,i+mat.KU+1); j++){
00093 operator()(i,j)-=mat(i,j);
00094 }
00095 }
00096
00097 mat.destroy();
00098 return *this;
00099 }
00100 else{
00101 zgbmatrix newmat(M,N,max(KL,mat.KL),max(KU,mat.KU));
00102 newmat.zero();
00103 for(long i=0; i<M; i++){
00104 for(long j=max(0,i-KL); j<min(N,i+KU+1); j++){
00105 newmat(i,j)+=operator()(i,j);
00106 }
00107 for(long j=max(0,i-mat.KL); j<min(mat.N,i+mat.KU+1); j++){
00108 newmat(i,j)-=mat(i,j);
00109 }
00110 }
00111
00112 swap(*this,newmat);
00113 mat.destroy();
00114 return *this;
00115 }
00116 }
00117
00118
00120 inline zgbmatrix& zgbmatrix::operator*=(const _zgbmatrix& mat)
00121 {
00122 #ifdef CPPL_VERBOSE
00123 std::cerr << "# [MARK] zgbmatrix::operator*=(const _zgbmatrix&)"
00124 << std::endl;
00125 #endif//CPPL_VERBOSE
00126
00127 #ifdef CPPL_DEBUG
00128 if(N!=mat.M){
00129 std::cerr << "[ERROR] zgbmatrix::operator*=(_zgbmatrix&)" << std::endl
00130 << "These two matrises can not make a product." << std::endl
00131 << "Your input was (" << M << "x" << N << ") * ("
00132 << mat.M << "x" << mat.N << ")." << std::endl;
00133 exit(1);
00134 }
00135 #endif//CPPL_DEBUG
00136
00137 zgbmatrix newmat( M, mat.N, KU+mat.KL-1, KL+mat.KU+1 );
00138 newmat.zero();
00139
00140 for(long i=0; i<newmat.M; i++){
00141 for(long j=max(0,i-newmat.KL); j<min(newmat.N,i+newmat.KU+1); j++){
00142 for(long k=max( max(0,i-KL), max(0,j-mat.KU) );
00143 k< min( min(N,i+KU+1), min(mat.M,j+mat.KL+1) ); k++){
00144 newmat(i,j)+= operator()(i,k)*mat(k,j);
00145 }
00146 }
00147 }
00148
00149 swap(*this,newmat);
00150 mat.destroy();
00151 return *this;
00152 }
00153
00157
00158
00160 inline _zgbmatrix operator+(const zgbmatrix& matA, const _zgbmatrix& matB)
00161 {
00162 #ifdef CPPL_VERBOSE
00163 std::cerr << "# [MARK] operator+(const zgbmatrix&, const _zgbmatrix&)"
00164 << std::endl;
00165 #endif//CPPL_VERBOSE
00166
00167 #ifdef CPPL_DEBUG
00168 if(matA.N!=matB.N || matA.M!=matB.M){
00169 std::cerr << "[ERROR] operator+(zgbmatrix&, _zgbmatrix&)" << std::endl
00170 << "These two matrises can not make a summation." << std::endl
00171 << "Your input was (" << matA.M << "x" << matA.N << ") + ("
00172 << matB.M << "x" << matB.N << ")." << std::endl;
00173 exit(1);
00174 }
00175 #endif//CPPL_DEBUG
00176
00177 if(matB.KL>matA.KL && matB.KU>matA.KU){
00178 for(long i=0; i<matA.M; i++){
00179 for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
00180 matB(i,j)+=matA(i,j);
00181 }
00182 }
00183
00184 return matB;
00185 }
00186 else{
00187 zgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
00188 newmat.zero();
00189
00190 for(long i=0; i<matA.M; i++){
00191 for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
00192 newmat(i,j)+=matA(i,j);
00193 }
00194 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00195 newmat(i,j)+=matB(i,j);
00196 }
00197 }
00198
00199 matB.destroy();
00200 return _(newmat);
00201 }
00202 }
00203
00204
00206 inline _zgbmatrix operator-(const zgbmatrix& matA, const _zgbmatrix& matB)
00207 {
00208 #ifdef CPPL_VERBOSE
00209 std::cerr << "# [MARK] operator-(const zgbmatrix&, const _zgbmatrix&)"
00210 << std::endl;
00211 #endif//CPPL_VERBOSE
00212
00213 #ifdef CPPL_DEBUG
00214 if(matA.N!=matB.N || matA.M!=matB.M){
00215 std::cerr << "[ERROR] operator-(zgbmatrix&, _zgbmatrix&)" << std::endl
00216 << "These two matrises can not make a subtraction." << std::endl
00217 << "Your input was (" << matA.M << "x" << matA.N << ") - ("
00218 << matB.M << "x" << matB.N << ")." << std::endl;
00219 exit(1);
00220 }
00221 #endif//CPPL_DEBUG
00222
00223 zgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
00224 newmat.zero();
00225
00226 for(long i=0; i<matA.M; i++){
00227 for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
00228 newmat(i,j)+=matA(i,j);
00229 }
00230 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00231 newmat(i,j)-=matB(i,j);
00232 }
00233 }
00234
00235 matB.destroy();
00236 return _(newmat);
00237 }
00238
00239
00241 inline _zgbmatrix operator*(const zgbmatrix& matA, const _zgbmatrix& matB)
00242 {
00243 #ifdef CPPL_VERBOSE
00244 std::cerr << "# [MARK] operator*(const zgbmatrix&, const _zgbmatrix&)"
00245 << std::endl;
00246 #endif//CPPL_VERBOSE
00247
00248 #ifdef CPPL_DEBUG
00249 if(matA.N!=matB.M){
00250 std::cerr << "[ERROR] operator*(zgbmatrix&, _zgbmatrix&)" << std::endl
00251 << "These two matrises can not make a product." << std::endl
00252 << "Your input was (" << matA.M << "x" << matA.N << ") * ("
00253 << matB.M << "x" << matB.N << ")." << std::endl;
00254 exit(1);
00255 }
00256 #endif//CPPL_DEBUG
00257
00258 zgbmatrix newmat( matA.M, matB.N, matA.KU+matB.KL-1, matA.KL+matB.KU+1 );
00259 newmat.zero();
00260
00261 long i, j, k;
00262 #pragma omp parallel for private(j,k)
00263 for(i=0; i<newmat.m; i++){
00264 for(j=max(0,i-newmat.kl); j<min(newmat.n,i+newmat.ku+1); j++){
00265 for(k=max( max(0,i-matA.KL), max(0,j-matB.KU) );
00266 k< min( min(matA.N,i+matA.KU+1), min(matB.M,j+matB.KL+1) ); k++){
00267 newmat(i,j)+= matA(i,k)*matB(k,j);
00268 }
00269 }
00270 }
00271
00272 matB.destroy();
00273 return _(newmat);
00274 }