00001
00003 inline _dgematrix operator+(const _dsymatrix& matA, const dgbmatrix& matB)
00004 {
00005 #ifdef CPPL_VERBOSE
00006 std::cerr << "# [MARK] operator+(const _dsymatrix&, const dgbmatrix&)"
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+(_dgematrix&, dgbmatrix&)" << 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 dgematrix newmat(matA.N, matA.N);
00021 for(long i=0; i<matB.M; i++){
00022 for(long j=i; j<matA.N; j++){
00023 newmat(i,j) = newmat(j,i) = 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 return _(newmat);
00032 }
00033
00034
00036 inline _dgematrix operator-(const _dsymatrix& matA, const dgbmatrix& matB)
00037 {
00038 #ifdef CPPL_VERBOSE
00039 std::cerr << "# [MARK] operator-(const _dsymatrix&, const dgbmatrix&)"
00040 << std::endl;
00041 #endif//CPPL_VERBOSE
00042
00043 #ifdef CPPL_DEBUG
00044 if(matA.N!=matB.N || matA.N!=matB.M){
00045 std::cerr << "[ERROR] operator+(_dsymatrix&, dgbmatrix&)" << std::endl
00046 << "These two matrises can not make a summation." << std::endl
00047 << "Your input was (" << matA.N << "x" << matA.N << ") + ("
00048 << matB.M << "x" << matB.N << ")." << std::endl;
00049 exit(1);
00050 }
00051 #endif//CPPL_DEBUG
00052
00053 dgematrix newmat(matA.N, matA.N);
00054 for(long i=0; i<matB.M; i++){
00055 for(long j=i; j<matA.N; j++){
00056 newmat(i,j) = newmat(j,i) = matA(i,j);
00057 }
00058 for(long j=max(0,i-matB.KL); j<min(matB.N,i+matB.KU+1); j++){
00059 newmat(i,j)-=matB(i,j);
00060 }
00061 }
00062
00063 matA.destroy();
00064 return _(newmat);
00065 }
00066
00067
00069 inline _dgematrix operator*(const _dsymatrix& matA, const dgbmatrix& matB)
00070 {
00071 #ifdef CPPL_VERBOSE
00072 std::cerr << "# [MARK] operator*(const _dsymatrix&, const dgbmatrix&)"
00073 << std::endl;
00074 #endif//CPPL_VERBOSE
00075
00076 #ifdef CPPL_DEBUG
00077 if(matA.N!=matB.M){
00078 std::cerr << "[ERROR] operator*(_dsymatrix&, dgbmatrix&)" << std::endl
00079 << "These two matrises can not make a product." << std::endl
00080 << "Your input was (" << matA.N << "x" << matA.N << ") * ("
00081 << matB.M << "x" << matB.N << ")." << std::endl;
00082 exit(1);
00083 }
00084 #endif//CPPL_DEBUG
00085
00086 dgematrix newmat( matA.N, matB.N );
00087 newmat.zero();
00088
00089 long i, j, k;
00090 #pragma omp parallel for private(j,k)
00091 for(i=0; i<newmat.m; i++){
00092 for(j=0; j<newmat.n; j++){
00093 for(k=max(0,j-matB.KU); k<min(matB.M,j+matB.KL+1); k++){
00094 newmat(i,j)+=matA(i,k)*matB(k,j);
00095 }
00096 }
00097 }
00098
00099 matA.destroy();
00100 return _(newmat);
00101 }