- Generated by
1.7.5.1
|
C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
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
1.7.5.1