C-XSC - A C++ Class Library for Extended Scientific Computing: imatrix.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: imatrix.cpp,v 1.25 2011/06/07 15:17:39 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