C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
rmatrix.cpp
00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
00007 **                          Universitaet Wuppertal, Germany   
00008 **
00009 **  This library is free software; you can redistribute it and/or
00010 **  modify it under the terms of the GNU Library General Public
00011 **  License as published by the Free Software Foundation; either
00012 **  version 2 of the License, or (at your option) any later version.
00013 **
00014 **  This library is distributed in the hope that it will be useful,
00015 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
00016 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 **  Library General Public License for more details.
00018 **
00019 **  You should have received a copy of the GNU Library General Public
00020 **  License along with this library; if not, write to the Free
00021 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 /* CVS $Id: rmatrix.cpp,v 1.28 2014/01/30 17:23:48 cxsc Exp $ */
00025 
00026 #define _CXSC_CPP
00027 
00028 #include "rmatrix.hpp"
00029 #include "matrix.inl"
00030 #include "rmatrix.inl"
00031 
00032 #include "dotk.inl"
00033 
00034 namespace cxsc {
00035   
00036     //Ostrowskis comparison matrix
00037     rmatrix CompMat ( const rmatrix& A) throw() {
00038       rmatrix M(Lb(A,1), Ub(A,1), Lb(A,2), Ub(A,2));
00039 
00040       for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) {
00041         for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) {
00042           if(i-Lb(A,1) == j-Lb(A,2))
00043             M[i][j] = abs(A[i][j]);
00044           else
00045             M[i][j] = -abs(A[i][j]);
00046         }
00047       }
00048 
00049   return M;
00050 }
00051 
00052       rmatrix Id ( const rmatrix& A )                        // Real identity matrix
00053       {                                                      //---------------------
00054         int i,j;
00055         int lbi = Lb(A,1), ubi = Ub(A,1);
00056         int lbj = Lb(A,2), ubj = Ub(A,2);
00057         rmatrix B(lbi,ubi,lbj,ubj);
00058       
00059         for (i = lbi; i <= ubi; i++)
00060           for (j = lbj; j <= ubj; j++)
00061             B[i][j] = (i==j) ? 1.0 : 0.0;
00062         return B;
00063       }
00064       
00065       rmatrix transp ( const rmatrix& A )                       // Transposed matrix
00066       {                                                         //------------------
00067         int      n;
00068         rmatrix  res(Lb(A,2),Ub(A,2),Lb(A,1),Ub(A,1));
00069       
00070         for (n = Lb(A,1); n <= Ub(A,1); n++) Col(res,n) = Row(A,n);
00071         return res;
00072       }
00073 
00074       void DoubleSize ( rmatrix& A )
00075       {
00076         int n = Lb(A,1);
00077         Resize(A,n,2*Ub(A,1)-n+1,Lb(A,2),Ub(A,2));
00078       }
00079 
00080 
00081         void accumulate(dotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00082 #if(CXSC_INDEX_CHECK)
00083         throw(OP_WITH_WRONG_DIM)
00084 #else
00085         throw()
00086 #endif
00087         { 
00088 #if(CXSC_INDEX_CHECK)
00089                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rmatrix_subv &, const rmatrix_subv &)"));
00090 #endif
00091                 addDot(dp,rv1,rv2);
00092         }
00093 
00094         void accumulate_approx(dotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2) {
00095                 addDot_op(dp,rv1,rv2);
00096         }
00097 
00098 
00099          void accumulate(dotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
00100 #if(CXSC_INDEX_CHECK)
00101         throw(OP_WITH_WRONG_DIM)
00102 #else
00103         throw()
00104 #endif
00105         { 
00106 #if(CXSC_INDEX_CHECK)
00107                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rvector &, const rmatrix_subv &)"));
00108 #endif
00109                 addDot(dp,rv1,rv2);
00110         }
00111 
00112         void accumulate_approx(dotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2) {
00113                 addDot_op(dp,rv1,rv2);
00114         }
00115 
00116 
00117         void accumulate(dotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
00118 #if(CXSC_INDEX_CHECK)
00119         throw(OP_WITH_WRONG_DIM)
00120 #else
00121         throw()
00122 #endif
00123         { 
00124 #if(CXSC_INDEX_CHECK)
00125                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rmatrix_subv &, const rvector &)"));
00126 #endif
00127                 addDot(dp,rv1,rv2);
00128         }
00129 
00130         void accumulate_approx(dotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2) {
00131                 addDot_op(dp,rv1,rv2);
00132         }
00133 
00134 
00135          void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00136 #if(CXSC_INDEX_CHECK)
00137         throw(OP_WITH_WRONG_DIM)
00138 #else
00139         throw()
00140 #endif
00141         { 
00142 #if(CXSC_INDEX_CHECK)
00143                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rmatrix_subv &, const rmatrix_subv &)"));
00144 #endif
00145                 dotprecision tmp(0.0);
00146                 tmp.set_k(dp.get_k());
00147                 addDot(tmp,rv1,rv2);
00148                 dp += tmp;
00149         }
00150 
00151 
00152 
00153          void accumulate(idotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
00154 #if(CXSC_INDEX_CHECK)
00155         throw(OP_WITH_WRONG_DIM)
00156 #else
00157         throw()
00158 #endif
00159         { 
00160 #if(CXSC_INDEX_CHECK)
00161                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rvector&, const rmatrix_subv &)"));
00162 #endif
00163                 dotprecision tmp(0.0);
00164                 tmp.set_k(dp.get_k());
00165                 addDot(tmp,rv1,rv2);
00166                 dp += tmp;
00167         }
00168 
00169 
00170 
00171          void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
00172 #if(CXSC_INDEX_CHECK)
00173         throw(OP_WITH_WRONG_DIM)
00174 #else
00175         throw()
00176 #endif
00177         { 
00178 #if(CXSC_INDEX_CHECK)
00179                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rmatrix_subv&, const rvector &)"));
00180 #endif
00181                 dotprecision tmp(0.0);
00182                 tmp.set_k(dp.get_k());
00183                 addDot(tmp,rv1,rv2);
00184                 dp += tmp;
00185         }
00186 
00187 
00188 
00189          void accumulate(cdotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00190 #if(CXSC_INDEX_CHECK)
00191         throw(OP_WITH_WRONG_DIM)
00192 #else
00193         throw()
00194 #endif
00195         { 
00196 #if(CXSC_INDEX_CHECK)
00197                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rmatrix_subv&, const rmatrix_subv &)"));
00198 #endif
00199                 addDot(Re(dp),rv1,rv2);
00200         }
00201 
00202          void accumulate_approx(cdotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00203         { 
00204                 addDot_op(Re(dp),rv1,rv2);
00205         }
00206 
00207 
00208          void accumulate(cdotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
00209 #if(CXSC_INDEX_CHECK)
00210         throw(OP_WITH_WRONG_DIM)
00211 #else
00212         throw()
00213 #endif
00214         { 
00215 #if(CXSC_INDEX_CHECK)
00216                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rvector&, const rmatrix_subv &)"));
00217 #endif
00218                 addDot(Re(dp),rv1,rv2);
00219         }
00220 
00221          void accumulate_approx(cdotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
00222         { 
00223                 addDot_op(Re(dp),rv1,rv2);
00224         }
00225 
00226          void accumulate(cdotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
00227 #if(CXSC_INDEX_CHECK)
00228         throw(OP_WITH_WRONG_DIM)
00229 #else
00230         throw()
00231 #endif
00232         { 
00233 #if(CXSC_INDEX_CHECK)
00234                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rmatrix_subv&, const rvector &)"));
00235 #endif
00236                 addDot(Re(dp),rv1,rv2);
00237         }
00238 
00239          void accumulate_approx(cdotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
00240         { 
00241                 addDot_op(Re(dp),rv1,rv2);
00242         }
00243 
00244          void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00245 #if(CXSC_INDEX_CHECK)
00246         throw(OP_WITH_WRONG_DIM)
00247 #else
00248         throw()
00249 #endif
00250         { 
00251 #if(CXSC_INDEX_CHECK)
00252                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const rmatrix_subv &)"));
00253 #endif
00254                 dotprecision tmp(0.0);
00255                 tmp.set_k(dp.get_k());
00256                 addDot(tmp,rv1,rv2);
00257                 dp += tmp;
00258         }
00259 
00260          void accumulate(cidotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
00261 #if(CXSC_INDEX_CHECK)
00262         throw(OP_WITH_WRONG_DIM)
00263 #else
00264         throw()
00265 #endif
00266         { 
00267 #if(CXSC_INDEX_CHECK)
00268                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rvector &, const rmatrix_subv &)"));
00269 #endif
00270                 dotprecision tmp(0.0);
00271                 tmp.set_k(dp.get_k());
00272                 addDot(tmp,rv1,rv2);
00273                 dp += tmp;
00274         }
00275 
00276          void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
00277 #if(CXSC_INDEX_CHECK)
00278         throw(OP_WITH_WRONG_DIM)
00279 #else
00280         throw()
00281 #endif
00282         { 
00283 #if(CXSC_INDEX_CHECK)
00284                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const rvector &)"));
00285 #endif
00286                 dotprecision tmp(0.0);
00287                 tmp.set_k(dp.get_k());
00288                 addDot(tmp,rv1,rv2);
00289                 dp += tmp;
00290         }
00291 
00292 
00293         void accumulate(dotprecision &dp,const rvector_slice &sl,const rmatrix_subv &sv)
00294 #if(CXSC_INDEX_CHECK)
00295         throw(OP_WITH_WRONG_DIM)
00296 #else
00297         throw()
00298 #endif
00299         {
00300 #if(CXSC_INDEX_CHECK)
00301                 if(VecLen(sl)!=VecLen(sv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rvector_slice &, const rmatrix_subv &)"));
00302 #endif
00303                 addDot(dp,sl,sv);
00304         }
00305 
00306         void accumulate_approx(dotprecision &dp,const rvector_slice &sl,const rmatrix_subv &sv) {
00307                 addDot_op(dp,sl,sv);
00308         }
00309 
00310 
00311          void accumulate(cdotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
00312 #if(CXSC_INDEX_CHECK)
00313         throw(OP_WITH_WRONG_DIM)
00314 #else
00315         throw()
00316 #endif
00317         { 
00318 #if(CXSC_INDEX_CHECK)
00319                 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rvector_slice&, const rmatrix_subv &)"));
00320 #endif  
00321                 addDot(Re(dp),sl1,rv2);
00322         }
00323 
00324          void accumulate_approx(cdotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
00325         { 
00326                 addDot_op(Re(dp),sl1,rv2);
00327         }
00328 
00329 
00330          void accumulate(idotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
00331 #if(CXSC_INDEX_CHECK)
00332         throw(OP_WITH_WRONG_DIM)
00333 #else
00334         throw()
00335 #endif
00336         { 
00337 #if(CXSC_INDEX_CHECK)
00338                 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rvector_slice&, const rmatrix_subv &)"));
00339 #endif
00340                 dotprecision tmp(0.0);
00341                 tmp.set_k(dp.get_k());
00342                 addDot(tmp,sl1,rv2);
00343                 dp += tmp;
00344         }
00345 
00346 
00347          void accumulate(cidotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
00348 #if(CXSC_INDEX_CHECK)
00349         throw(OP_WITH_WRONG_DIM)
00350 #else
00351         throw()
00352 #endif
00353         { 
00354 #if(CXSC_INDEX_CHECK)
00355                 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rvector_slice &, const rmatrix_subv &)"));
00356 #endif
00357                 dotprecision tmp(0.0);
00358                 tmp.set_k(dp.get_k());
00359                 addDot(tmp,sl1,rv2);
00360                 dp += tmp;
00361         }
00362 
00363 
00364         void accumulate(dotprecision &dp,const rmatrix_subv &mv,const rvector_slice &vs)
00365 #if(CXSC_INDEX_CHECK)
00366         throw(OP_WITH_WRONG_DIM)
00367 #else
00368         throw()
00369 #endif
00370         { 
00371 #if(CXSC_INDEX_CHECK)
00372                 if(VecLen(mv)!=VecLen(vs)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rmatrix_subv &, const rvector_slice &)"));
00373 #endif
00374                 addDot(dp,vs,mv); 
00375         }
00376 
00377         void accumulate_approx(dotprecision &dp,const rmatrix_subv &mv,const rvector_slice &vs) {
00378                 addDot_op(dp,vs,mv);
00379         }
00380 
00381 
00382 
00383          void accumulate(cdotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
00384 #if(CXSC_INDEX_CHECK)
00385         throw(OP_WITH_WRONG_DIM)
00386 #else
00387         throw()
00388 #endif
00389         { 
00390 #if(CXSC_INDEX_CHECK)
00391                 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rmatrix_subv&, const rvector_slice &)"));
00392 #endif  
00393                 addDot(Re(dp),rv1,sl2);
00394         }
00395 
00396          void accumulate_approx(cdotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
00397         { 
00398                 addDot_op(Re(dp),rv1,sl2);
00399         }
00400 
00401 
00402          void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
00403 #if(CXSC_INDEX_CHECK)
00404         throw(OP_WITH_WRONG_DIM)
00405 #else
00406         throw()
00407 #endif
00408         { 
00409 #if(CXSC_INDEX_CHECK)
00410                 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rmatrix_subv&, const rvector_slice &)"));
00411 #endif
00412                 dotprecision tmp(0.0);
00413                 tmp.set_k(dp.get_k());
00414                 addDot(tmp,rv1,sl2);
00415                 dp += tmp;
00416         }
00417 
00418 
00419          void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
00420 #if(CXSC_INDEX_CHECK)
00421         throw(OP_WITH_WRONG_DIM)
00422 #else
00423         throw()
00424 #endif
00425         { 
00426 #if(CXSC_INDEX_CHECK)
00427                 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const rvector_slice &)"));
00428 #endif
00429                 dotprecision tmp(0.0);
00430                 tmp.set_k(dp.get_k());
00431                 addDot(tmp,rv1,sl2);
00432                 dp += tmp;
00433         }
00434 
00435 
00436 } // namespace cxsc