00001
00003 inline _zgematrix operator+(const _zhematrix& matA, const _zgbmatrix& matB)
00004 {
00005 #ifdef CPPL_VERBOSE
00006 std::cerr << "# [MARK] operator+(const _zhematrix&, const _zgbmatrix&)"
00007 << std::endl;
00008 #endif//CPPL_VERBOSE
00009
00010 #ifdef CPPL_DEBUG
00011 if(matA.N!=matB.N || matA.N!=matB.M){
00012 std::cerr << "[ERROR] operator+(_zhematrix&, _zgbmatrix&)" << std::endl
00013 << "These two matrises can not make a summation." << std::endl
00014 << "Your input was (" << matA.N << "x" << matA.N << ") + ("
00015 << matB.M << "x" << matB.N << ")." << std::endl;
00016 exit(1);
00017 }
00018 #endif//CPPL_DEBUG
00019
00020 zgematrix newmat(matA.N, matA.N);
00021 for(long i=0; i<matB.M; i++){
00022 for(long j=0; j<matA.N; j++){
00023 newmat(i,j) =matA(i,j);
00024 }
00025 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00026 newmat(i,j)+=matB(i,j);
00027 }
00028 }
00029
00030 matA.destroy();
00031 matB.destroy();
00032 return _(newmat);
00033 }
00034
00035
00037 inline _zgematrix operator-(const _zhematrix& matA, const _zgbmatrix& matB)
00038 {
00039 #ifdef CPPL_VERBOSE
00040 std::cerr << "# [MARK] operator-(const _zhematrix&, const _zgbmatrix&)"
00041 << std::endl;
00042 #endif//CPPL_VERBOSE
00043
00044 #ifdef CPPL_DEBUG
00045 if(matA.N!=matB.N || matA.N!=matB.M){
00046 std::cerr << "[ERROR] operator+(_zhematrix&, _zgbmatrix&)" << std::endl
00047 << "These two matrises can not make a summation." << std::endl
00048 << "Your input was (" << matA.N << "x" << matA.N << ") + ("
00049 << matB.M << "x" << matB.N << ")." << std::endl;
00050 exit(1);
00051 }
00052 #endif//CPPL_DEBUG
00053
00054 zgematrix newmat(matA.N, matA.N);
00055 for(long i=0; i<matB.M; i++){
00056 for(long j=0; j<matA.N; j++){
00057 newmat(i,j) =matA(i,j);
00058 }
00059 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00060 newmat(i,j)-=matB(i,j);
00061 }
00062 }
00063
00064 matA.destroy();
00065 matB.destroy();
00066 return _(newmat);
00067 }
00068
00069
00071 inline _zgematrix operator*(const _zhematrix& matA, const _zgbmatrix& matB)
00072 {
00073 #ifdef CPPL_VERBOSE
00074 std::cerr << "# [MARK] operator*(const _zhematrix&, const _zgbmatrix&)"
00075 << std::endl;
00076 #endif//CPPL_VERBOSE
00077
00078 #ifdef CPPL_DEBUG
00079 if(matA.N!=matB.M){
00080 std::cerr << "[ERROR] operator*(_zhematrix&, _zgbmatrix&)" << std::endl
00081 << "These two matrises can not make a product." << std::endl
00082 << "Your input was (" << matA.N << "x" << matA.N << ") * ("
00083 << matB.M << "x" << matB.N << ")." << std::endl;
00084 exit(1);
00085 }
00086 #endif//CPPL_DEBUG
00087
00088 zgematrix newmat( matA.N, matB.N );
00089 newmat.zero();
00090
00091 long i, j, k;
00092 #pragma omp parallel for private(j,k)
00093 for(i=0; i<newmat.m; i++){
00094 for(j=0; j<newmat.n; j++){
00095 for(k=max(0,j-matB.KU); k<min(matB.M,j+matB.KL+1); k++){
00096 newmat(i,j)+=matA(i,k)*matB(k,j);
00097 }
00098 }
00099 }
00100
00101 matA.destroy();
00102 matB.destroy();
00103 return _(newmat);
00104 }