00001
00003 inline _dgbmatrix operator+(const _dgbmatrix& matA, const _dgbmatrix& matB)
00004 {
00005 #ifdef CPPL_VERBOSE
00006 std::cerr << "# [MARK] operator+(const _dgbmatrix&, const _dgbmatrix&)"
00007 << std::endl;
00008 #endif//CPPL_VERBOSE
00009
00010 #ifdef CPPL_DEBUG
00011 if(matA.N!=matB.N || matA.M!=matB.M){
00012 std::cerr << "[ERROR] operator+(_dgbmatrix&, _dgbmatrix&)" << std::endl
00013 << "These two matrises can not make a summation." << std::endl
00014 << "Your input was (" << matA.M << "x" << matA.N << ") + ("
00015 << matB.M << "x" << matB.N << ")." << std::endl;
00016 exit(1);
00017 }
00018 #endif//CPPL_DEBUG
00019
00020 if(matA.KL>matB.KL && matA.KU>matB.KU){
00021 for(long i=0; i<matB.M; i++){
00022 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00023 matA(i,j)+=matB(i,j);
00024 }
00025 }
00026
00027 matB.destroy();
00028 return matA;
00029 }
00030
00031 else if(matB.KL>matA.KL && matB.KU>matA.KU){
00032 for(long i=0; i<matA.M; i++){
00033 for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
00034 matB(i,j)+=matA(i,j);
00035 }
00036 }
00037
00038 matA.destroy();
00039 return matB;
00040 }
00041
00042 else{
00043 dgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
00044 newmat.zero();
00045
00046 for(long i=0; i<matA.M; i++){
00047 for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
00048 newmat(i,j)+=matA(i,j);
00049 }
00050 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00051 newmat(i,j)+=matB(i,j);
00052 }
00053 }
00054
00055 matA.destroy();
00056 matB.destroy();
00057 return _(newmat);
00058 }
00059 }
00060
00061
00063 inline _dgbmatrix operator-(const _dgbmatrix& matA, const _dgbmatrix& matB)
00064 {
00065 #ifdef CPPL_VERBOSE
00066 std::cerr << "# [MARK] operator-(const _dgbmatrix&, const _dgbmatrix&)"
00067 << std::endl;
00068 #endif//CPPL_VERBOSE
00069
00070 #ifdef CPPL_DEBUG
00071 if(matA.N!=matB.N || matA.M!=matB.M){
00072 std::cerr << "[ERROR] operator-(_dgbmatrix&, _dgbmatrix&)" << std::endl
00073 << "These two matrises can not make a subtraction." << std::endl
00074 << "Your input was (" << matA.M << "x" << matA.N << ") - ("
00075 << matB.M << "x" << matB.N << ")." << std::endl;
00076 exit(1);
00077 }
00078 #endif//CPPL_DEBUG
00079
00080 if(matA.KL>matB.KL && matA.KU>matB.KU){
00081 for(long i=0; i<matB.M; i++){
00082 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00083 matA(i,j)-=matB(i,j);
00084 }
00085 }
00086
00087 matB.destroy();
00088 return matA;
00089 }
00090
00091 else{
00092 dgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
00093 newmat.zero();
00094
00095 for(long i=0; i<matA.M; i++){
00096 for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
00097 newmat(i,j)+=matA(i,j);
00098 }
00099 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00100 newmat(i,j)-=matB(i,j);
00101 }
00102 }
00103
00104 matA.destroy();
00105 matB.destroy();
00106 return _(newmat);
00107 }
00108 }
00109
00110
00112 inline _dgbmatrix operator*(const _dgbmatrix& matA, const _dgbmatrix& matB)
00113 {
00114 #ifdef CPPL_VERBOSE
00115 std::cerr << "# [MARK] operator*(const _dgbmatrix&, const _dgbmatrix&)"
00116 << std::endl;
00117 #endif//CPPL_VERBOSE
00118
00119 #ifdef CPPL_DEBUG
00120 if(matA.N!=matB.M){
00121 std::cerr << "[ERROR] operator*(_dgbmatrix&, _dgbmatrix&)" << std::endl
00122 << "These two matrises can not make a product." << std::endl
00123 << "Your input was (" << matA.M << "x" << matA.N << ") * ("
00124 << matB.M << "x" << matB.N << ")." << std::endl;
00125 exit(1);
00126 }
00127 #endif//CPPL_DEBUG
00128
00129 dgbmatrix newmat( matA.M, matB.N, matA.KU+matB.KL-1, matA.KL+matB.KU+1 );
00130 newmat.zero();
00131
00132 long i, j, k;
00133 #pragma omp parallel for private(j,k)
00134 for(i=0; i<newmat.m; i++){
00135 for(j=max(0,i-newmat.kl); j<min(newmat.n,i+newmat.ku+1); j++){
00136 for(k=max( max(0,i-matA.KL), max(0,j-matB.KU) );
00137 k< min( min(matA.N,i+matA.KU+1), min(matB.M,j+matB.KL+1) ); k++){
00138 newmat(i,j)+= matA(i,k)*matB(k,j);
00139 }
00140 }
00141 }
00142
00143 matA.destroy();
00144 matB.destroy();
00145 return _(newmat);
00146 }