C-XSC - A C++ Class Library for Extended Scientific Computing: rmatrix.cpp Source File

00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.1)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2011 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.23 2011/06/07 15:17:41 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       rmatrix Id ( const rmatrix& A )                        // Real identity matrix
00037       {                                                      //---------------------
00038         int i,j;
00039         int lbi = Lb(A,1), ubi = Ub(A,1);
00040         int lbj = Lb(A,2), ubj = Ub(A,2);
00041         rmatrix B(lbi,ubi,lbj,ubj);
00042       
00043         for (i = lbi; i <= ubi; i++)
00044           for (j = lbj; j <= ubj; j++)
00045             B[i][j] = (i==j) ? 1.0 : 0.0;
00046         return B;
00047       }
00048       
00049       rmatrix transp ( const rmatrix& A )                       // Transposed matrix
00050       {                                                         //------------------
00051         int      n;
00052         rmatrix  res(Lb(A,2),Ub(A,2),Lb(A,1),Ub(A,1));
00053       
00054         for (n = Lb(A,1); n <= Ub(A,1); n++) Col(res,n) = Row(A,n);
00055         return res;
00056       }
00057 
00058       void DoubleSize ( rmatrix& A )
00059       {
00060         int n = Lb(A,1);
00061         Resize(A,n,2*Ub(A,1)-n+1,Lb(A,2),Ub(A,2));
00062       }
00063 
00064 
00065         void accumulate(dotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00066 #if(CXSC_INDEX_CHECK)
00067         throw(OP_WITH_WRONG_DIM)
00068 #else
00069         throw()
00070 #endif
00071         { 
00072 #if(CXSC_INDEX_CHECK)
00073                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rmatrix_subv &, const rmatrix_subv &)"));
00074 #endif
00075                 addDot(dp,rv1,rv2);
00076         }
00077 
00078         void accumulate_approx(dotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2) {
00079                 addDot_op(dp,rv1,rv2);
00080         }
00081 
00082 
00083          void accumulate(dotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
00084 #if(CXSC_INDEX_CHECK)
00085         throw(OP_WITH_WRONG_DIM)
00086 #else
00087         throw()
00088 #endif
00089         { 
00090 #if(CXSC_INDEX_CHECK)
00091                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rvector &, const rmatrix_subv &)"));
00092 #endif
00093                 addDot(dp,rv1,rv2);
00094         }
00095 
00096         void accumulate_approx(dotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2) {
00097                 addDot_op(dp,rv1,rv2);
00098         }
00099 
00100 
00101         void accumulate(dotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
00102 #if(CXSC_INDEX_CHECK)
00103         throw(OP_WITH_WRONG_DIM)
00104 #else
00105         throw()
00106 #endif
00107         { 
00108 #if(CXSC_INDEX_CHECK)
00109                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rmatrix_subv &, const rvector &)"));
00110 #endif
00111                 addDot(dp,rv1,rv2);
00112         }
00113 
00114         void accumulate_approx(dotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2) {
00115                 addDot_op(dp,rv1,rv2);
00116         }
00117 
00118 
00119          void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00120 #if(CXSC_INDEX_CHECK)
00121         throw(OP_WITH_WRONG_DIM)
00122 #else
00123         throw()
00124 #endif
00125         { 
00126 #if(CXSC_INDEX_CHECK)
00127                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rmatrix_subv &, const rmatrix_subv &)"));
00128 #endif
00129                 dotprecision tmp(0.0);
00130                 tmp.set_k(dp.get_k());
00131                 addDot(tmp,rv1,rv2);
00132                 dp += tmp;
00133         }
00134 
00135 
00136 
00137          void accumulate(idotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
00138 #if(CXSC_INDEX_CHECK)
00139         throw(OP_WITH_WRONG_DIM)
00140 #else
00141         throw()
00142 #endif
00143         { 
00144 #if(CXSC_INDEX_CHECK)
00145                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rvector&, const rmatrix_subv &)"));
00146 #endif
00147                 dotprecision tmp(0.0);
00148                 tmp.set_k(dp.get_k());
00149                 addDot(tmp,rv1,rv2);
00150                 dp += tmp;
00151         }
00152 
00153 
00154 
00155          void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
00156 #if(CXSC_INDEX_CHECK)
00157         throw(OP_WITH_WRONG_DIM)
00158 #else
00159         throw()
00160 #endif
00161         { 
00162 #if(CXSC_INDEX_CHECK)
00163                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rmatrix_subv&, const rvector &)"));
00164 #endif
00165                 dotprecision tmp(0.0);
00166                 tmp.set_k(dp.get_k());
00167                 addDot(tmp,rv1,rv2);
00168                 dp += tmp;
00169         }
00170 
00171 
00172 
00173          void accumulate(cdotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00174 #if(CXSC_INDEX_CHECK)
00175         throw(OP_WITH_WRONG_DIM)
00176 #else
00177         throw()
00178 #endif
00179         { 
00180 #if(CXSC_INDEX_CHECK)
00181                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rmatrix_subv&, const rmatrix_subv &)"));
00182 #endif
00183                 addDot(Re(dp),rv1,rv2);
00184         }
00185 
00186          void accumulate_approx(cdotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00187         { 
00188                 addDot_op(Re(dp),rv1,rv2);
00189         }
00190 
00191 
00192          void accumulate(cdotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
00193 #if(CXSC_INDEX_CHECK)
00194         throw(OP_WITH_WRONG_DIM)
00195 #else
00196         throw()
00197 #endif
00198         { 
00199 #if(CXSC_INDEX_CHECK)
00200                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rvector&, const rmatrix_subv &)"));
00201 #endif
00202                 addDot(Re(dp),rv1,rv2);
00203         }
00204 
00205          void accumulate_approx(cdotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
00206         { 
00207                 addDot_op(Re(dp),rv1,rv2);
00208         }
00209 
00210          void accumulate(cdotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
00211 #if(CXSC_INDEX_CHECK)
00212         throw(OP_WITH_WRONG_DIM)
00213 #else
00214         throw()
00215 #endif
00216         { 
00217 #if(CXSC_INDEX_CHECK)
00218                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rmatrix_subv&, const rvector &)"));
00219 #endif
00220                 addDot(Re(dp),rv1,rv2);
00221         }
00222 
00223          void accumulate_approx(cdotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
00224         { 
00225                 addDot_op(Re(dp),rv1,rv2);
00226         }
00227 
00228          void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
00229 #if(CXSC_INDEX_CHECK)
00230         throw(OP_WITH_WRONG_DIM)
00231 #else
00232         throw()
00233 #endif
00234         { 
00235 #if(CXSC_INDEX_CHECK)
00236                 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const rmatrix_subv &)"));
00237 #endif
00238                 dotprecision tmp(0.0);
00239                 tmp.set_k(dp.get_k());
00240                 addDot(tmp,rv1,rv2);
00241                 dp += tmp;
00242         }
00243 
00244          void accumulate(cidotprecision &dp, const rvector & 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 rvector &, 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 rmatrix_subv & rv1, const rvector &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 rmatrix_subv &, const rvector &)"));
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 
00277         void accumulate(dotprecision &dp,const rvector_slice &sl,const rmatrix_subv &sv)
00278 #if(CXSC_INDEX_CHECK)
00279         throw(OP_WITH_WRONG_DIM)
00280 #else
00281         throw()
00282 #endif
00283         {
00284 #if(CXSC_INDEX_CHECK)
00285                 if(VecLen(sl)!=VecLen(sv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rvector_slice &, const rmatrix_subv &)"));
00286 #endif
00287                 addDot(dp,sl,sv);
00288         }
00289 
00290         void accumulate_approx(dotprecision &dp,const rvector_slice &sl,const rmatrix_subv &sv) {
00291                 addDot_op(dp,sl,sv);
00292         }
00293 
00294 
00295          void accumulate(cdotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
00296 #if(CXSC_INDEX_CHECK)
00297         throw(OP_WITH_WRONG_DIM)
00298 #else
00299         throw()
00300 #endif
00301         { 
00302 #if(CXSC_INDEX_CHECK)
00303                 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rvector_slice&, const rmatrix_subv &)"));
00304 #endif  
00305                 addDot(Re(dp),sl1,rv2);
00306         }
00307 
00308          void accumulate_approx(cdotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
00309         { 
00310                 addDot_op(Re(dp),sl1,rv2);
00311         }
00312 
00313 
00314          void accumulate(idotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
00315 #if(CXSC_INDEX_CHECK)
00316         throw(OP_WITH_WRONG_DIM)
00317 #else
00318         throw()
00319 #endif
00320         { 
00321 #if(CXSC_INDEX_CHECK)
00322                 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rvector_slice&, const rmatrix_subv &)"));
00323 #endif
00324                 dotprecision tmp(0.0);
00325                 tmp.set_k(dp.get_k());
00326                 addDot(tmp,sl1,rv2);
00327                 dp += tmp;
00328         }
00329 
00330 
00331          void accumulate(cidotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
00332 #if(CXSC_INDEX_CHECK)
00333         throw(OP_WITH_WRONG_DIM)
00334 #else
00335         throw()
00336 #endif
00337         { 
00338 #if(CXSC_INDEX_CHECK)
00339                 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rvector_slice &, const rmatrix_subv &)"));
00340 #endif
00341                 dotprecision tmp(0.0);
00342                 tmp.set_k(dp.get_k());
00343                 addDot(tmp,sl1,rv2);
00344                 dp += tmp;
00345         }
00346 
00347 
00348         void accumulate(dotprecision &dp,const rmatrix_subv &mv,const rvector_slice &vs)
00349 #if(CXSC_INDEX_CHECK)
00350         throw(OP_WITH_WRONG_DIM)
00351 #else
00352         throw()
00353 #endif
00354         { 
00355 #if(CXSC_INDEX_CHECK)
00356                 if(VecLen(mv)!=VecLen(vs)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rmatrix_subv &, const rvector_slice &)"));
00357 #endif
00358                 addDot(dp,vs,mv); 
00359         }
00360 
00361         void accumulate_approx(dotprecision &dp,const rmatrix_subv &mv,const rvector_slice &vs) {
00362                 addDot_op(dp,vs,mv);
00363         }
00364 
00365 
00366 
00367          void accumulate(cdotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
00368 #if(CXSC_INDEX_CHECK)
00369         throw(OP_WITH_WRONG_DIM)
00370 #else
00371         throw()
00372 #endif
00373         { 
00374 #if(CXSC_INDEX_CHECK)
00375                 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rmatrix_subv&, const rvector_slice &)"));
00376 #endif  
00377                 addDot(Re(dp),rv1,sl2);
00378         }
00379 
00380          void accumulate_approx(cdotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
00381         { 
00382                 addDot_op(Re(dp),rv1,sl2);
00383         }
00384 
00385 
00386          void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
00387 #if(CXSC_INDEX_CHECK)
00388         throw(OP_WITH_WRONG_DIM)
00389 #else
00390         throw()
00391 #endif
00392         { 
00393 #if(CXSC_INDEX_CHECK)
00394                 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rmatrix_subv&, const rvector_slice &)"));
00395 #endif
00396                 dotprecision tmp(0.0);
00397                 tmp.set_k(dp.get_k());
00398                 addDot(tmp,rv1,sl2);
00399                 dp += tmp;
00400         }
00401 
00402 
00403          void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
00404 #if(CXSC_INDEX_CHECK)
00405         throw(OP_WITH_WRONG_DIM)
00406 #else
00407         throw()
00408 #endif
00409         { 
00410 #if(CXSC_INDEX_CHECK)
00411                 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const rvector_slice &)"));
00412 #endif
00413                 dotprecision tmp(0.0);
00414                 tmp.set_k(dp.get_k());
00415                 addDot(tmp,rv1,sl2);
00416                 dp += tmp;
00417         }
00418 
00419 
00420 } // namespace cxsc