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: ivector.cpp,v 1.26 2014/01/30 17:23:45 cxsc Exp $ */ 00025 00026 #define _CXSC_CPP 00027 00028 #include "ivector.hpp" 00029 #include "vector.inl" 00030 #include "ivector.inl" 00031 00032 #include "idotk.inl" 00033 00034 namespace cxsc { 00035 00036 int in ( const ivector& x, const ivector& y ) // Inner inclusion for two ivectors 00037 { //--------------------------------- 00038 int i = Lb(x), n = Ub(x), ok = 1; 00039 00040 while (ok && i <= n) { ok = in(x[i],y[i]); i++; } 00041 return ok; 00042 } 00043 00044 int in ( int x, ivector& y ) // Inner inclusion of an integer in a ivector 00045 { //------------------------------------------- 00046 int i = Lb(y), n = Ub(y), ok = 1; 00047 double d = x; 00048 00049 while (ok && i <= n) { ok = in(d,y[i]); i++; } 00050 return ok; 00051 } 00052 00053 ivector Blow ( const ivector& x, real eps ) // Epsilon deflation 00054 { //------------------ 00055 int i; 00056 ivector h(Lb(x),Ub(x)); 00057 00058 for (i = Lb(x); i <= Ub(x); i++) h[i] = Blow(x[i],eps); 00059 return h; 00060 } 00061 00062 int Disjoint ( ivector& a, ivector& b ) // Test for disjointedness 00063 { //------------------------ 00064 int al = Lb(a), au = Ub(a), bl = Lb(b); 00065 int disjointed, i, d; 00066 00067 d = bl - al; 00068 00069 i = al; 00070 disjointed = 0; 00071 do { 00072 if (Disjoint(a[i],b[i+d])) 00073 disjointed = 1; 00074 else 00075 i++; 00076 } while ( !(disjointed || i > au) ); 00077 return disjointed; 00078 } 00079 00080 int Zero ( ivector& x ) // Check for zero vector 00081 { //---------------------- 00082 int i, ok; 00083 for (i = Lb(x), ok = 1; i <= Ub(x) && ok; i++) ok = (x[i] == 0.0); 00084 return ok; 00085 } 00086 00087 rvector mid ( ivector& v ) // Vector of midpoints 00088 { //-------------------- 00089 int i; 00090 int l = Lb(v), u = Ub(v); 00091 rvector w(l,u); 00092 00093 for (i = l; i <= u; i++) w[i] = mid(v[i]); 00094 return w; 00095 } 00096 00097 real MaxRelDiam ( const ivector& v ) // Maximum relative diameter 00098 { //-------------------------- 00099 real r; 00100 int i, l=Lb(v), u=Ub(v); 00101 00102 r = 0.0; 00103 for (i = l; i <= u; i++) 00104 if (RelDiam(v[i]) > r) r = RelDiam(v[i]); 00105 return r; 00106 } 00107 00108 real MaxRelDiam ( const ivector_slice& v ) // Maximum relative diameter 00109 { //-------------------------- 00110 real r; 00111 int i, l=Lb(v), u=Ub(v); 00112 00113 r = 0.0; 00114 for (i = l; i <= u; i++) 00115 if (RelDiam(v[i]) > r) r = RelDiam(v[i]); 00116 return r; 00117 } 00118 00119 //---------------------------------------------------------------------------- 00120 // Checks if the diameter of the interval vector 'v' is less or equal to 'n' 00121 // ulps. Ulp is an abbreviation for: units in the last place of the mantissa. 00122 //---------------------------------------------------------------------------- 00123 int UlpAcc ( ivector& v, int n ) 00124 { 00125 int k, upper; 00126 00127 for (k = Lb(v), upper = Ub(v); (k < upper) && UlpAcc(v[k],n); k++); 00128 return UlpAcc(v[k],n); 00129 } 00130 00131 // The 'DoubleSize' functions double the number of rows of a matrix 00132 // or double the length of a vector preserving existing components. 00133 //------------------------------------------------------------------ 00134 void DoubleSize ( ivector& x ) 00135 { 00136 int n = Lb(x); 00137 Resize(x,n,2*Ub(x)-n+1); 00138 } 00139 00140 00141 //accurate dot product definitions 00142 00143 00144 void accumulate(idotprecision &dp, const ivector & rv1, const ivector &rv2) 00145 #if(CXSC_INDEX_CHECK) 00146 throw(OP_WITH_WRONG_DIM) 00147 #else 00148 throw() 00149 #endif 00150 { 00151 #if(CXSC_INDEX_CHECK) 00152 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, ivector &)")); 00153 #endif 00154 addDot(dp,rv1,rv2); 00155 } 00156 00157 00158 void accumulate(idotprecision &dp, const ivector_slice & sl, const ivector &rv) 00159 #if(CXSC_INDEX_CHECK) 00160 throw(OP_WITH_WRONG_DIM) 00161 #else 00162 throw() 00163 #endif 00164 { 00165 #if(CXSC_INDEX_CHECK) 00166 if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, ivector &)")); 00167 #endif 00168 addDot(dp,sl,rv); 00169 } 00170 00171 00172 void accumulate(idotprecision &dp, const ivector &rv, const ivector_slice &sl) 00173 #if(CXSC_INDEX_CHECK) 00174 throw(OP_WITH_WRONG_DIM) 00175 #else 00176 throw() 00177 #endif 00178 { 00179 #if(CXSC_INDEX_CHECK) 00180 if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, ivector_slice &)")); 00181 #endif 00182 addDot(dp,rv,sl); 00183 } 00184 00185 00186 void accumulate(idotprecision &dp, const ivector_slice & sl1, const ivector_slice &sl2) 00187 #if(CXSC_INDEX_CHECK) 00188 throw(OP_WITH_WRONG_DIM) 00189 #else 00190 throw() 00191 #endif 00192 { 00193 #if(CXSC_INDEX_CHECK) 00194 if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, ivector_slice &)")); 00195 #endif 00196 addDot(dp,sl1,sl2); 00197 } 00198 00199 00200 void accumulate(cidotprecision &dp, const ivector & rv1, const ivector &rv2) 00201 #if(CXSC_INDEX_CHECK) 00202 throw(OP_WITH_WRONG_DIM) 00203 #else 00204 throw() 00205 #endif 00206 { 00207 #if(CXSC_INDEX_CHECK) 00208 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector &, ivector &)")); 00209 #endif 00210 idotprecision tmp = Re(dp); 00211 tmp.set_k(dp.get_k()); 00212 addDot(tmp,rv1,rv2); 00213 SetRe(dp,tmp); 00214 } 00215 00216 00217 void accumulate(cidotprecision &dp, const ivector_slice & sl, const ivector &rv) 00218 #if(CXSC_INDEX_CHECK) 00219 throw(OP_WITH_WRONG_DIM) 00220 #else 00221 throw() 00222 #endif 00223 { 00224 #if(CXSC_INDEX_CHECK) 00225 if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector_slice &, ivector &)")); 00226 #endif 00227 idotprecision tmp = Re(dp); 00228 tmp.set_k(dp.get_k()); 00229 addDot(tmp,sl,rv); 00230 SetRe(dp,tmp); 00231 } 00232 00233 00234 void accumulate(cidotprecision &dp, const ivector &rv, const ivector_slice &sl) 00235 #if(CXSC_INDEX_CHECK) 00236 throw(OP_WITH_WRONG_DIM) 00237 #else 00238 throw() 00239 #endif 00240 { 00241 #if(CXSC_INDEX_CHECK) 00242 if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector &, ivector_slice &)")); 00243 #endif 00244 idotprecision tmp = Re(dp); 00245 tmp.set_k(dp.get_k()); 00246 addDot(tmp,rv,sl); 00247 SetRe(dp,tmp); 00248 } 00249 00250 00251 void accumulate(cidotprecision &dp, const ivector_slice & sl1, const ivector_slice &sl2) 00252 #if(CXSC_INDEX_CHECK) 00253 throw(OP_WITH_WRONG_DIM) 00254 #else 00255 throw() 00256 #endif 00257 { 00258 #if(CXSC_INDEX_CHECK) 00259 if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector_slice &, ivector_slice &)")); 00260 #endif 00261 idotprecision tmp = Re(dp); 00262 tmp.set_k(dp.get_k()); 00263 addDot(tmp,sl1,sl2); 00264 SetRe(dp,tmp); 00265 } 00266 00267 00268 void accumulate(idotprecision &dp, const rvector & rv1, const ivector &rv2) 00269 #if(CXSC_INDEX_CHECK) 00270 throw(OP_WITH_WRONG_DIM) 00271 #else 00272 throw() 00273 #endif 00274 { 00275 #if(CXSC_INDEX_CHECK) 00276 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector &, ivector &)")); 00277 #endif 00278 addDot(dp,rv1,rv2); 00279 } 00280 00281 00282 void accumulate(idotprecision &dp, const ivector & rv1, const rvector &rv2) 00283 #if(CXSC_INDEX_CHECK) 00284 throw(OP_WITH_WRONG_DIM) 00285 #else 00286 throw() 00287 #endif 00288 { 00289 #if(CXSC_INDEX_CHECK) 00290 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, rvector &)")); 00291 #endif 00292 addDot(dp,rv1,rv2); 00293 } 00294 00295 00296 void accumulate(idotprecision &dp, const rvector_slice & sl, const ivector &rv) 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(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector_slice &, ivector &)")); 00305 #endif 00306 addDot(dp,sl,rv); 00307 } 00308 00309 00310 void accumulate(idotprecision &dp,const ivector_slice &sl,const rvector &rv) 00311 #if(CXSC_INDEX_CHECK) 00312 throw(OP_WITH_WRONG_DIM) 00313 #else 00314 throw() 00315 #endif 00316 { 00317 #if(CXSC_INDEX_CHECK) 00318 if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, rvector &)")); 00319 #endif 00320 addDot(dp,sl,rv); 00321 } 00322 00323 00324 void accumulate(idotprecision &dp, const rvector &rv, const ivector_slice &sl) 00325 #if(CXSC_INDEX_CHECK) 00326 throw(OP_WITH_WRONG_DIM) 00327 #else 00328 throw() 00329 #endif 00330 { 00331 #if(CXSC_INDEX_CHECK) 00332 if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector &, ivector_slice &)")); 00333 #endif 00334 addDot(dp,rv,sl); 00335 } 00336 00337 00338 void accumulate(idotprecision &dp,const ivector &rv,const rvector_slice &sl) 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(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, rvector_slice &)")); 00347 #endif 00348 addDot(dp,rv,sl); 00349 } 00350 00351 00352 void accumulate(idotprecision &dp, const ivector_slice & sl1, const rvector_slice &sl2) 00353 #if(CXSC_INDEX_CHECK) 00354 throw(OP_WITH_WRONG_DIM) 00355 #else 00356 throw() 00357 #endif 00358 { 00359 #if(CXSC_INDEX_CHECK) 00360 if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, rvector_slice &)")); 00361 #endif 00362 addDot(dp,sl1,sl2); 00363 } 00364 00365 00366 void accumulate(idotprecision &dp, const rvector_slice & sl1, const ivector_slice &sl2) 00367 #if(CXSC_INDEX_CHECK) 00368 throw(OP_WITH_WRONG_DIM) 00369 #else 00370 throw() 00371 #endif 00372 { 00373 #if(CXSC_INDEX_CHECK) 00374 if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector_slice &, ivector_slice &)")); 00375 #endif 00376 addDot(dp,sl1,sl2); 00377 } 00378 00379 00380 void accumulate(cidotprecision &dp, const rvector & rv1, const ivector &rv2) 00381 #if(CXSC_INDEX_CHECK) 00382 throw(OP_WITH_WRONG_DIM) 00383 #else 00384 throw() 00385 #endif 00386 { 00387 #if(CXSC_INDEX_CHECK) 00388 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const rvector &, ivector &)")); 00389 #endif 00390 idotprecision tmp = Re(dp); 00391 tmp.set_k(dp.get_k()); 00392 addDot(tmp,rv1,rv2); 00393 SetRe(dp,tmp); 00394 } 00395 00396 00397 void accumulate(cidotprecision &dp, const ivector & rv1, const rvector &rv2) 00398 #if(CXSC_INDEX_CHECK) 00399 throw(OP_WITH_WRONG_DIM) 00400 #else 00401 throw() 00402 #endif 00403 { 00404 #if(CXSC_INDEX_CHECK) 00405 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector &, rvector &)")); 00406 #endif 00407 idotprecision tmp = Re(dp); 00408 tmp.set_k(dp.get_k()); 00409 addDot(tmp,rv1,rv2); 00410 SetRe(dp,tmp); 00411 } 00412 00413 00414 void accumulate(cidotprecision &dp, const rvector_slice & sl, const ivector &rv) 00415 #if(CXSC_INDEX_CHECK) 00416 throw(OP_WITH_WRONG_DIM) 00417 #else 00418 throw() 00419 #endif 00420 { 00421 #if(CXSC_INDEX_CHECK) 00422 if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const rvector_slice &, ivector &)")); 00423 #endif 00424 idotprecision tmp = Re(dp); 00425 tmp.set_k(dp.get_k()); 00426 addDot(tmp,sl,rv); 00427 SetRe(dp,tmp); 00428 } 00429 00430 00431 void accumulate(cidotprecision &dp,const ivector_slice &sl,const rvector &rv) 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(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector_slice &, rvector &)")); 00440 #endif 00441 idotprecision tmp = Re(dp); 00442 tmp.set_k(dp.get_k()); 00443 addDot(tmp,sl,rv); 00444 SetRe(dp,tmp); 00445 } 00446 00447 00448 void accumulate(cidotprecision &dp, const rvector &rv, const ivector_slice &sl) 00449 #if(CXSC_INDEX_CHECK) 00450 throw(OP_WITH_WRONG_DIM) 00451 #else 00452 throw() 00453 #endif 00454 { 00455 #if(CXSC_INDEX_CHECK) 00456 if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const rvector &, ivector_slice &)")); 00457 #endif 00458 idotprecision tmp = Re(dp); 00459 tmp.set_k(dp.get_k()); 00460 addDot(tmp,rv,sl); 00461 SetRe(dp,tmp); 00462 } 00463 00464 00465 void accumulate(cidotprecision &dp,const ivector &rv,const rvector_slice &sl) 00466 #if(CXSC_INDEX_CHECK) 00467 throw(OP_WITH_WRONG_DIM) 00468 #else 00469 throw() 00470 #endif 00471 { 00472 #if(CXSC_INDEX_CHECK) 00473 if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector &, rvector_slice &)")); 00474 #endif 00475 idotprecision tmp = Re(dp); 00476 tmp.set_k(dp.get_k()); 00477 addDot(tmp,rv,sl); 00478 SetRe(dp,tmp); 00479 } 00480 00481 00482 void accumulate(cidotprecision &dp, const ivector_slice & sl1, const rvector_slice &sl2) 00483 #if(CXSC_INDEX_CHECK) 00484 throw(OP_WITH_WRONG_DIM) 00485 #else 00486 throw() 00487 #endif 00488 { 00489 #if(CXSC_INDEX_CHECK) 00490 if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector_slice &, rvector_slice &)")); 00491 #endif 00492 idotprecision tmp = Re(dp); 00493 tmp.set_k(dp.get_k()); 00494 addDot(tmp,sl1,sl2); 00495 SetRe(dp,tmp); 00496 } 00497 00498 00499 void accumulate(cidotprecision &dp, const rvector_slice & sl1, const ivector_slice &sl2) 00500 #if(CXSC_INDEX_CHECK) 00501 throw(OP_WITH_WRONG_DIM) 00502 #else 00503 throw() 00504 #endif 00505 { 00506 #if(CXSC_INDEX_CHECK) 00507 if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const rvector_slice &, ivector_slice &)")); 00508 #endif 00509 idotprecision tmp = Re(dp); 00510 tmp.set_k(dp.get_k()); 00511 addDot(tmp,sl1,sl2); 00512 SetRe(dp,tmp); 00513 } 00514 00515 00516 //Summation accumulates 00517 void accumulate(idotprecision &dp, const ivector& v) { 00518 addSum(Inf(dp),Inf(v)); 00519 addSum(Sup(dp),Sup(v)); 00520 } 00521 00522 void accumulate(cidotprecision &dp, const ivector& v) { 00523 addSum(InfRe(dp),Inf(v)); 00524 addSum(SupRe(dp),Sup(v)); 00525 } 00526 00527 } // namespace cxsc 00528