- 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: imatrix.cpp,v 1.29 2014/01/30 17:23:45 cxsc Exp $ */ 00025 00026 #define _CXSC_CPP 00027 00028 #include "imatrix.hpp" 00029 #include "matrix.inl" 00030 #include "imatrix.inl" 00031 #include "ivecrmat.inl" 00032 00033 #include "idotk.inl" 00034 00035 namespace cxsc { 00036 00037 //Ostrowskis comparison matrix 00038 rmatrix CompMat ( const imatrix& A) { 00039 rmatrix M(Lb(A,1), Ub(A,1), Lb(A,2), Ub(A,2)); 00040 00041 for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) { 00042 for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) { 00043 if(i-Lb(A,1) == j-Lb(A,2)) 00044 M[i][j] = AbsMin(A[i][j]); 00045 else 00046 M[i][j] = -AbsMax(A[i][j]); 00047 } 00048 } 00049 00050 return M; 00051 } 00052 00053 00054 imatrix Id ( const imatrix& A ) // Interval identity matrix 00055 { //------------------------- 00056 int i,j; 00057 int lbi = Lb(A,1), ubi = Ub(A,1); 00058 int lbj = Lb(A,2), ubj = Ub(A,2); 00059 imatrix B(lbi,ubi,lbj,ubj); 00060 00061 for (i = lbi; i <= ubi; i++) 00062 for (j = lbj; j <= ubj; j++) 00063 B[i][j] = interval( (i==j) ? 1.0 : 0.0 ); 00064 return B; 00065 } 00066 00067 imatrix transp ( const imatrix& A ) // Transposed matrix 00068 { //------------------ 00069 int n; 00070 imatrix res(Lb(A,2),Ub(A,2),Lb(A,1),Ub(A,1)); 00071 00072 for (n = Lb(A,1); n <= Ub(A,1); n++) Col(res,n) = Row(A,n); 00073 return res; 00074 } 00075 00076 real MaxRelDiam ( const imatrix_subv& v ) // Maximum relative diameter 00077 { //-------------------------- 00078 real r; 00079 int i, l=Lb(v), u=Ub(v); 00080 00081 r = 0.0; 00082 for (i = l; i <= u; i++) 00083 if (RelDiam(v[i]) > r) r = RelDiam(v[i]); 00084 return r; 00085 } 00086 00087 // The 'DoubleSize' functions double the number of rows of a matrix 00088 // or double the length of a vector preserving existing components. 00089 //------------------------------------------------------------------ 00090 00091 void DoubleSize ( imatrix& A ) 00092 { 00093 int n = Lb(A,1); 00094 Resize(A,n,2*Ub(A,1)-n+1,Lb(A,2),Ub(A,2)); 00095 } 00096 00097 00098 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const imatrix_subv &rv2) 00099 #if(CXSC_INDEX_CHECK) 00100 throw(OP_WITH_WRONG_DIM) 00101 #else 00102 throw() 00103 #endif 00104 { 00105 #if(CXSC_INDEX_CHECK) 00106 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, imatrix_subv &)")); 00107 #endif 00108 addDot(dp,rv1,rv2); 00109 } 00110 00111 void accumulate(idotprecision &dp, const ivector & rv1, const imatrix_subv &rv2) 00112 #if(CXSC_INDEX_CHECK) 00113 throw(OP_WITH_WRONG_DIM) 00114 #else 00115 throw() 00116 #endif 00117 { 00118 #if(CXSC_INDEX_CHECK) 00119 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, imatrix_subv &)")); 00120 #endif 00121 addDot(dp,rv1,rv2); 00122 } 00123 00124 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const ivector &rv2) 00125 #if(CXSC_INDEX_CHECK) 00126 throw(OP_WITH_WRONG_DIM) 00127 #else 00128 throw() 00129 #endif 00130 { 00131 #if(CXSC_INDEX_CHECK) 00132 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, ivector &)")); 00133 #endif 00134 addDot(dp,rv1,rv2); 00135 } 00136 00137 00138 void accumulate(idotprecision &dp, const ivector_slice & sl1, const imatrix_subv &rv2) 00139 #if(CXSC_INDEX_CHECK) 00140 throw(OP_WITH_WRONG_DIM) 00141 #else 00142 throw() 00143 #endif 00144 { 00145 #if(CXSC_INDEX_CHECK) 00146 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, imatrix_subv &)")); 00147 #endif 00148 addDot(dp,sl1,rv2); 00149 } 00150 00151 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const ivector_slice &sl2) 00152 #if(CXSC_INDEX_CHECK) 00153 throw(OP_WITH_WRONG_DIM) 00154 #else 00155 throw() 00156 #endif 00157 { 00158 #if(CXSC_INDEX_CHECK) 00159 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, ivector_slice &)")); 00160 #endif 00161 addDot(dp,rv1,sl2); 00162 } 00163 00164 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const imatrix_subv &rv2) 00165 #if(CXSC_INDEX_CHECK) 00166 throw(OP_WITH_WRONG_DIM) 00167 #else 00168 throw() 00169 #endif 00170 { 00171 #if(CXSC_INDEX_CHECK) 00172 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const imatrix_subv &)")); 00173 #endif 00174 idotprecision tmp(0.0); 00175 tmp.set_k(dp.get_k()); 00176 addDot(tmp,rv1,rv2); 00177 dp += tmp; 00178 } 00179 00180 void accumulate(cidotprecision &dp, const ivector & rv1, const imatrix_subv &rv2) 00181 #if(CXSC_INDEX_CHECK) 00182 throw(OP_WITH_WRONG_DIM) 00183 #else 00184 throw() 00185 #endif 00186 { 00187 #if(CXSC_INDEX_CHECK) 00188 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const ivector &, const imatrix_subv &)")); 00189 #endif 00190 idotprecision tmp(0.0); 00191 tmp.set_k(dp.get_k()); 00192 addDot(tmp,rv1,rv2); 00193 dp += tmp; 00194 } 00195 00196 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const ivector &rv2) 00197 #if(CXSC_INDEX_CHECK) 00198 throw(OP_WITH_WRONG_DIM) 00199 #else 00200 throw() 00201 #endif 00202 { 00203 #if(CXSC_INDEX_CHECK) 00204 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const ivector &)")); 00205 #endif 00206 idotprecision tmp(0.0); 00207 tmp.set_k(dp.get_k()); 00208 addDot(tmp,rv1,rv2); 00209 dp += tmp; 00210 } 00211 00212 void accumulate(cidotprecision &dp, const ivector_slice & sl1, const imatrix_subv &rv2) 00213 #if(CXSC_INDEX_CHECK) 00214 throw(OP_WITH_WRONG_DIM) 00215 #else 00216 throw() 00217 #endif 00218 { 00219 #if(CXSC_INDEX_CHECK) 00220 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const ivector_slice &, const imatrix_subv &)")); 00221 #endif 00222 idotprecision tmp(0.0); 00223 tmp.set_k(dp.get_k()); 00224 addDot(tmp,sl1,rv2); 00225 dp += tmp; 00226 } 00227 00228 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const ivector_slice &sl2) 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(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const ivector_slice &)")); 00237 #endif 00238 idotprecision tmp(0.0); 00239 tmp.set_k(dp.get_k()); 00240 addDot(tmp,rv1,sl2); 00241 dp += tmp; 00242 } 00243 00244 void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const imatrix_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(idotprecision &, const rmatrix_subv &, imatrix_subv &)")); 00253 #endif 00254 addDot(dp,rv1,rv2); 00255 } 00256 00257 void accumulate(idotprecision &dp, const rvector & rv1, const imatrix_subv &rv2) 00258 #if(CXSC_INDEX_CHECK) 00259 throw(OP_WITH_WRONG_DIM) 00260 #else 00261 throw() 00262 #endif 00263 { 00264 #if(CXSC_INDEX_CHECK) 00265 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector &, imatrix_subv &)")); 00266 #endif 00267 addDot(dp,rv1,rv2); 00268 } 00269 00270 void accumulate(idotprecision &dp, const rvector_slice & sl1, const imatrix_subv &rv2) 00271 #if(CXSC_INDEX_CHECK) 00272 throw(OP_WITH_WRONG_DIM) 00273 #else 00274 throw() 00275 #endif 00276 { 00277 #if(CXSC_INDEX_CHECK) 00278 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector_slice &, imatrix_subv &)")); 00279 #endif 00280 addDot(dp,sl1,rv2); 00281 } 00282 00283 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const rmatrix_subv &rv2) 00284 #if(CXSC_INDEX_CHECK) 00285 throw(OP_WITH_WRONG_DIM) 00286 #else 00287 throw() 00288 #endif 00289 { 00290 #if(CXSC_INDEX_CHECK) 00291 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, rmatrix_subv &)")); 00292 #endif 00293 addDot(dp,rv1,rv2); 00294 } 00295 00296 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const rvector &rv2) 00297 #if(CXSC_INDEX_CHECK) 00298 throw(OP_WITH_WRONG_DIM) 00299 #else 00300 throw() 00301 #endif 00302 { 00303 #if(CXSC_INDEX_CHECK) 00304 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, rvector &)")); 00305 #endif 00306 addDot(dp,rv1,rv2); 00307 } 00308 00309 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const rvector_slice &sl2) 00310 #if(CXSC_INDEX_CHECK) 00311 throw(OP_WITH_WRONG_DIM) 00312 #else 00313 throw() 00314 #endif 00315 { 00316 #if(CXSC_INDEX_CHECK) 00317 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, rvector_slice &)")); 00318 #endif 00319 addDot(dp,rv1,sl2); 00320 } 00321 00322 void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const imatrix_subv &rv2) 00323 #if(CXSC_INDEX_CHECK) 00324 throw(OP_WITH_WRONG_DIM) 00325 #else 00326 throw() 00327 #endif 00328 { 00329 #if(CXSC_INDEX_CHECK) 00330 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const imatrix_subv &)")); 00331 #endif 00332 idotprecision tmp(0.0); 00333 tmp.set_k(dp.get_k()); 00334 addDot(tmp,rv1,rv2); 00335 dp += tmp; 00336 } 00337 00338 void accumulate(cidotprecision &dp, const rvector & rv1, const imatrix_subv &rv2) 00339 #if(CXSC_INDEX_CHECK) 00340 throw(OP_WITH_WRONG_DIM) 00341 #else 00342 throw() 00343 #endif 00344 { 00345 #if(CXSC_INDEX_CHECK) 00346 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rvector &, const imatrix_subv &)")); 00347 #endif 00348 idotprecision tmp(0.0); 00349 tmp.set_k(dp.get_k()); 00350 addDot(tmp,rv1,rv2); 00351 dp += tmp; 00352 } 00353 00354 void accumulate(cidotprecision &dp, const rvector_slice & sl1, const imatrix_subv &rv2) 00355 #if(CXSC_INDEX_CHECK) 00356 throw(OP_WITH_WRONG_DIM) 00357 #else 00358 throw() 00359 #endif 00360 { 00361 #if(CXSC_INDEX_CHECK) 00362 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rvector_slice&, const imatrix_subv &)")); 00363 #endif 00364 idotprecision tmp(0.0); 00365 tmp.set_k(dp.get_k()); 00366 addDot(tmp,sl1,rv2); 00367 dp += tmp; 00368 } 00369 00370 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const rmatrix_subv &rv2) 00371 #if(CXSC_INDEX_CHECK) 00372 throw(OP_WITH_WRONG_DIM) 00373 #else 00374 throw() 00375 #endif 00376 { 00377 #if(CXSC_INDEX_CHECK) 00378 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const rmatrix_subv &)")); 00379 #endif 00380 idotprecision tmp(0.0); 00381 tmp.set_k(dp.get_k()); 00382 addDot(tmp,rv1,rv2); 00383 dp += tmp; 00384 } 00385 00386 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const rvector &rv2) 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(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const rvector &)")); 00395 #endif 00396 idotprecision tmp(0.0); 00397 tmp.set_k(dp.get_k()); 00398 addDot(tmp,rv1,rv2); 00399 dp += tmp; 00400 } 00401 00402 void accumulate(cidotprecision &dp, const imatrix_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(cidotprecision&, const imatrix_subv &, const rvector_slice &)")); 00411 #endif 00412 idotprecision tmp(0.0); 00413 tmp.set_k(dp.get_k()); 00414 addDot(tmp,rv1,sl2); 00415 dp += tmp; 00416 } 00417 00418 void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const ivector &rv2) 00419 #if(CXSC_INDEX_CHECK) 00420 throw(OP_WITH_WRONG_DIM) 00421 #else 00422 throw() 00423 #endif 00424 { 00425 #if(CXSC_INDEX_CHECK) 00426 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rmatrix_subv &, ivector &)")); 00427 #endif 00428 addDot(dp,rv1,rv2); 00429 } 00430 00431 void accumulate(idotprecision &dp, const ivector & rv1, const rmatrix_subv &rv2) 00432 #if(CXSC_INDEX_CHECK) 00433 throw(OP_WITH_WRONG_DIM) 00434 #else 00435 throw() 00436 #endif 00437 { 00438 #if(CXSC_INDEX_CHECK) 00439 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, rmatrix_subv &)")); 00440 #endif 00441 addDot(dp,rv1,rv2); 00442 } 00443 00444 void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const ivector &rv2) 00445 #if(CXSC_INDEX_CHECK) 00446 throw(OP_WITH_WRONG_DIM) 00447 #else 00448 throw() 00449 #endif 00450 { 00451 #if(CXSC_INDEX_CHECK) 00452 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const ivector &)")); 00453 #endif 00454 idotprecision tmp(0.0); 00455 tmp.set_k(dp.get_k()); 00456 addDot(tmp,rv1,rv2); 00457 dp += tmp; 00458 } 00459 00460 void accumulate(cidotprecision &dp, const ivector & rv1, const rmatrix_subv &rv2) 00461 #if(CXSC_INDEX_CHECK) 00462 throw(OP_WITH_WRONG_DIM) 00463 #else 00464 throw() 00465 #endif 00466 { 00467 #if(CXSC_INDEX_CHECK) 00468 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const ivector &, const rmatrix_subv &)")); 00469 #endif 00470 idotprecision tmp(0.0); 00471 tmp.set_k(dp.get_k()); 00472 addDot(tmp,rv1,rv2); 00473 dp += tmp; 00474 } 00475 00476 void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const ivector_slice &rv2) 00477 #if(CXSC_INDEX_CHECK) 00478 throw(OP_WITH_WRONG_DIM) 00479 #else 00480 throw() 00481 #endif 00482 { 00483 #if(CXSC_INDEX_CHECK) 00484 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rmatrix_subv &, ivector_slice &)")); 00485 #endif 00486 addDot(dp,rv1,rv2); 00487 } 00488 00489 void accumulate(idotprecision &dp, const ivector_slice & rv1, const rmatrix_subv &rv2) 00490 #if(CXSC_INDEX_CHECK) 00491 throw(OP_WITH_WRONG_DIM) 00492 #else 00493 throw() 00494 #endif 00495 { 00496 #if(CXSC_INDEX_CHECK) 00497 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, rmatrix_subv &)")); 00498 #endif 00499 addDot(dp,rv1,rv2); 00500 } 00501 00502 void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const ivector_slice &rv2) 00503 #if(CXSC_INDEX_CHECK) 00504 throw(OP_WITH_WRONG_DIM) 00505 #else 00506 throw() 00507 #endif 00508 { 00509 #if(CXSC_INDEX_CHECK) 00510 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const ivector_slice &)")); 00511 #endif 00512 idotprecision tmp(0.0); 00513 tmp.set_k(dp.get_k()); 00514 addDot(tmp,rv1,rv2); 00515 dp += tmp; 00516 } 00517 00518 void accumulate(cidotprecision &dp, const ivector_slice & rv1, const rmatrix_subv &rv2) 00519 #if(CXSC_INDEX_CHECK) 00520 throw(OP_WITH_WRONG_DIM) 00521 #else 00522 throw() 00523 #endif 00524 { 00525 #if(CXSC_INDEX_CHECK) 00526 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const ivector_slice &, const rmatrix_subv &)")); 00527 #endif 00528 idotprecision tmp(0.0); 00529 tmp.set_k(dp.get_k()); 00530 addDot(tmp,rv1,rv2); 00531 dp += tmp; 00532 } 00533 00534 00535 } // namespace cxsc 00536
1.7.5.1