00001
00003 inline _zgbmatrix operator+(const _zgbmatrix& matA, const zgbmatrix& matB)
00004 {
00005 #ifdef CPPL_VERBOSE
00006 std::cerr << "# [MARK] operator+(const _zgbmatrix&, const zgbmatrix&)"
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+(_zgbmatrix&, zgbmatrix&)" << 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 return matA;
00028 }
00029
00030 else{
00031 zgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
00032 newmat.zero();
00033
00034 for(long i=0; i<matA.M; i++){
00035 for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
00036 newmat(i,j)+=matA(i,j);
00037 }
00038 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00039 newmat(i,j)+=matB(i,j);
00040 }
00041 }
00042
00043 matA.destroy();
00044 return _(newmat);
00045 }
00046 }
00047
00048
00050 inline _zgbmatrix operator-(const _zgbmatrix& matA, const zgbmatrix& matB)
00051 {
00052 #ifdef CPPL_VERBOSE
00053 std::cerr << "# [MARK] operator-(const _zgbmatrix&, const zgbmatrix&)"
00054 << std::endl;
00055 #endif//CPPL_VERBOSE
00056
00057 #ifdef CPPL_DEBUG
00058 if(matA.N!=matB.N || matA.M!=matB.M){
00059 std::cerr << "[ERROR] operator-(_zgbmatrix&, zgbmatrix&)" << std::endl
00060 << "These two matrises can not make a subtraction." << std::endl
00061 << "Your input was (" << matA.M << "x" << matA.N << ") - ("
00062 << matB.M << "x" << matB.N << ")." << std::endl;
00063 exit(1);
00064 }
00065 #endif//CPPL_DEBUG
00066
00067 if(matA.KL>=matB.KL && matA.KU>=matB.KU){
00068 for(long i=0; i<matB.M; i++){
00069 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00070 matA(i,j)-=matB(i,j);
00071 }
00072 }
00073
00074 return matA;
00075 }
00076
00077 else{
00078 zgbmatrix newmat(matA.M,matA.N,max(matA.KL,matB.KL),max(matA.KU,matB.KU));
00079 newmat.zero();
00080
00081 for(long i=0; i<matA.M; i++){
00082 for(long j=max(0,i-matA.KL); j<min(matA.N,i+matA.KU+1); j++){
00083 newmat(i,j)+=matA(i,j);
00084 }
00085 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00086 newmat(i,j)-=matB(i,j);
00087 }
00088 }
00089
00090 return _(newmat);
00091 }
00092 }
00093
00094
00096 inline _zgbmatrix operator*(const _zgbmatrix& matA, const zgbmatrix& matB)
00097 {
00098 #ifdef CPPL_VERBOSE
00099 std::cerr << "# [MARK] operator*(const _zgbmatrix&, const zgbmatrix&)"
00100 << std::endl;
00101 #endif//CPPL_VERBOSE
00102
00103 #ifdef CPPL_DEBUG
00104 if(matA.N!=matB.M){
00105 std::cerr << "[ERROR] operator*(_zgbmatrix&, zgbmatrix&)" << std::endl
00106 << "These two matrises can not make a product." << std::endl
00107 << "Your input was (" << matA.M << "x" << matA.N << ") * ("
00108 << matB.M << "x" << matB.N << ")." << std::endl;
00109 exit(1);
00110 }
00111 #endif//CPPL_DEBUG
00112
00113 zgbmatrix newmat( matA.M, matB.N, matA.KU+matB.KL-1, matA.KL+matB.KU+1 );
00114 newmat.zero();
00115
00116 long i, j, k;
00117 #pragma omp parallel for private(j,k)
00118 for(i=0; i<newmat.m; i++){
00119 for(j=max(0,i-newmat.kl); j<min(newmat.n,i+newmat.ku+1); j++){
00120 for(k=max( max(0,i-matA.KL), max(0,j-matB.KU) );
00121 k< min( min(matA.N,i+matA.KU+1), min(matB.M,j+matB.KL+1) ); k++){
00122 newmat(i,j)+= matA(i,k)*matB(k,j);
00123 }
00124 }
00125 }
00126
00127 matA.destroy();
00128 return _(newmat);
00129 }