00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef INTERLIB_H
00018 #define INTERLIB_H
00019 #include <iostream>
00020
00021 using std::ostream;
00022 using std::istream;
00023
00030 namespace intblas
00031 {
00032
00035 enum NORMTYPE
00036 {
00037 NORM1,
00038 NORM2,
00039 NORMF,
00040 NORMMAX,
00041 NORMINF
00042 };
00043
00046 enum TRIMAT_TYPE {
00049 LTR_TOP,
00052 LTR_BOT,
00055 RTL_TOP,
00058 RTL_BOT
00059 };
00060
00061
00062 class IntervalMatrix;
00063 class IntervalVector;
00064
00110 struct INTERVAL_EXCEPTION
00111 {
00112 INTERVAL_EXCEPTION(){ msg[0] = 0; code = 0; }
00113
00114 char msg[200];
00115 int code;
00116 };
00117
00118 ostream& operator<<( ostream& out, const INTERVAL_EXCEPTION& ie );
00119
00139 class Interval
00140 {
00141 protected:
00142 double m_lower;
00143 double m_upper;
00144 double* m_ptr[2];
00145 static INTERVAL_EXCEPTION error;
00146
00147 public:
00148 static bool skip_roundout;
00149
00150 static void INIT_INTERVAL( bool skip_roundout = false );
00151
00152 public:
00153 Interval();
00154 Interval( double lower, double upper );
00155 Interval( double point );
00156 inline Interval( const Interval& i ){ this->operator=( i ); }
00157 ~Interval();
00158
00159
00160 static INTERVAL_EXCEPTION getLastError();
00161 static INTERVAL_EXCEPTION setError( int code, const char* msg = NULL );
00162
00163
00167 inline double& lowerBound(){ return m_lower; }
00171 inline double& upperBound(){ return m_upper; }
00176 inline double midpoint() const { return intblas::midpoint( *this ); }
00177
00182 inline Interval intervalHull( const Interval& b ) const
00183 { return intblas::intervalHull( *this, b ); }
00184
00189 inline double width() const { return intblas::width( *this ); }
00190
00191
00192 double& operator[]( int index ) const;
00193 friend ostream& operator<<( ostream& out, Interval const& a );
00194 friend istream& operator>>( istream& in, Interval& i );
00195 Interval& operator=( const Interval& b );
00196
00197
00198 Interval operator+( const Interval& b ) const;
00199 Interval operator-( const Interval& b ) const;
00200 Interval operator*( const Interval& b ) const;
00201 Interval operator/( const Interval& b ) const;
00202 Interval operator^( const Interval& y ) const;
00203 Interval operator*( double r ) const;
00204 friend Interval operator*( double r, const Interval& b ) { return (b * r); }
00205 Interval operator+( double r ) const;
00206 friend Interval operator+( double r, const Interval& b ) { return (b + r); }
00207 Interval operator-( double r ) const;
00208 friend Interval operator-( double r, const Interval& b );
00209 Interval operator/( double r ) const;
00210 friend Interval operator/( double r, const Interval& b );
00211 Interval operator-() const;
00212 Interval cancel( const Interval& b ) const;
00213 friend Interval cancel( const Interval &a, const Interval& b );
00214 Interval scalarAdd( double r ) const;
00215 Interval scalarMult( double r ) const;
00216
00217 void operator+=( const Interval& b );
00218 void operator+=( double r );
00219
00220 void operator-=( const Interval& b );
00221 void operator-=( double r );
00222
00223 void operator*=( const Interval& b );
00224 void operator*=( double r );
00225
00226 void operator/=( const Interval& b );
00227 void operator/=( double r );
00228
00229
00230 Interval operator&( const Interval& b ) const;
00231 Interval operator|( const Interval& b ) const;
00232 friend double midpoint( const Interval& a );
00233 friend Interval intervalHull( const Interval& a, const Interval& b );
00234 friend double width( const Interval& a );
00235
00236
00237 bool disjoint( const Interval& b ) const;
00238 bool encloses( const Interval& b ) const;
00239 bool encloses( double v ) const;
00240 bool interior( const Interval& b ) const;
00241 bool interior( double v ) const;
00242 bool operator==( const Interval& b ) const;
00243 bool operator==( const double v ) const;
00244 bool operator<( const Interval& b ) const;
00245 bool operator>( const Interval& b ) const;
00246 bool operator<( double b ) const;
00247 bool operator>( double b ) const;
00248
00249 bool isEmpty() const;
00250 friend bool isEmpty( const Interval& v );
00251
00252
00253 friend Interval isin( const Interval& x );
00254 friend Interval icos( const Interval& x );
00255 friend Interval itan( const Interval& x );
00256 friend Interval iasin( const Interval& xx );
00257 friend Interval iacos( const Interval& x );
00258 friend Interval iatan( const Interval& x );
00259
00260
00261 friend Interval isqrt( const Interval& x );
00262 friend Interval iexp( const Interval& x );
00263 friend Interval ilog( const Interval& x );
00264 friend Interval ipow( const Interval& x, int exp );
00265 friend Interval ipow( const Interval& x, const Interval& y );
00266 friend double iabs( const Interval& x );
00267 };
00268
00277 extern void roundOut( Interval& x, bool lower = true, bool upper = true );
00278
00285 class IntervalMatrix
00286 {
00287 protected:
00288 Interval **a;
00289 unsigned int m_cols;
00290 unsigned int m_rows;
00291 bool ir;
00292 Interval ***a_ptr;
00293 Interval *abase;
00294
00295 void alloc( unsigned int rows, unsigned int cols );
00296 public:
00297 static const Interval one;
00298
00299 public:
00303 inline IntervalMatrix(){ m_cols = 0; m_rows = 0; a = NULL; ir = false; };
00304 IntervalMatrix( unsigned int n, unsigned int m, bool is_ir = false );
00305 IntervalMatrix( const IntervalMatrix& m );
00306 IntervalMatrix( double** l, double** u, unsigned int rows, unsigned int cols );
00307 ~IntervalMatrix();
00308
00309
00313 inline int getRows() const { return (int)m_rows; }
00317 inline int getCols() const { return (int)m_cols; }
00325 inline Interval& getAt( unsigned int r, unsigned int c ) const { return a[r][c]; }
00332 inline void setAt( unsigned int r, unsigned int c, const Interval& i )
00333 { a[r][c] = i; }
00334
00335 IntervalMatrix getIdentity();
00336 friend IntervalMatrix getIdentity( unsigned int rows, unsigned int cols );
00337 static IntervalMatrix getIdentity( unsigned int rows, unsigned int cols );
00338
00339
00344 inline Interval* operator[]( unsigned int index ) const { return a[index]; };
00345 void operator=( const IntervalMatrix& b );
00346 bool operator==( const IntervalMatrix& b ) const;
00347
00348 friend void swap( IntervalMatrix& a, IntervalMatrix& b );
00349
00350 friend ostream& operator<<( ostream& out, const IntervalMatrix& mat );
00351 friend istream& operator>>( istream& in, IntervalMatrix& m );
00352
00353
00354 IntervalMatrix& operator+( IntervalMatrix& b );
00355 static void add( const IntervalMatrix& a, const IntervalMatrix& b,
00356 IntervalMatrix& result );
00357
00358 IntervalMatrix& operator+( const Interval& scalar );
00359 friend IntervalMatrix& operator+( const Interval& scalar, IntervalMatrix &a );
00360 static void add( const IntervalMatrix& a, const Interval& scalar,
00361 IntervalMatrix& result );
00362
00363 IntervalMatrix& operator-( IntervalMatrix& b );
00364 static void sub( const IntervalMatrix& a, const IntervalMatrix& b,
00365 IntervalMatrix& result );
00366
00367 IntervalMatrix& operator-( const Interval& scalar );
00368 friend IntervalMatrix& operator-( const Interval& scalar, IntervalMatrix &a );
00369 static void sub( IntervalMatrix& a, const Interval& scalar,
00370 IntervalMatrix& result, bool right = true );
00371
00372 IntervalMatrix& operator*( IntervalMatrix& b );
00373 static void mult( IntervalMatrix& a, IntervalMatrix& b,
00374 IntervalMatrix& result );
00375
00376 IntervalMatrix& operator*( const Interval& scalar );
00377 friend IntervalMatrix& operator*( const Interval& scalar, IntervalMatrix &a );
00378 static void mult( IntervalMatrix& a, const Interval& scalar,
00379 IntervalMatrix& result );
00380
00381 IntervalVector& operator*( IntervalVector& v );
00382 friend IntervalVector& operator*( IntervalVector& v, IntervalMatrix& a );
00383 static void mult( IntervalMatrix& a, IntervalVector& v,
00384 IntervalVector& result, bool right );
00385
00386
00387 IntervalMatrix& operator-();
00388 friend IntervalMatrix& neg( IntervalMatrix& a );
00389 static void neg( IntervalMatrix& a, IntervalMatrix& result );
00390
00391 IntervalMatrix& transpose();
00392 friend IntervalMatrix& transpose( IntervalMatrix& a );
00393 static void transpose( IntervalMatrix& a, IntervalMatrix& result );
00394
00395 double norm( NORMTYPE n );
00396 friend double norm( IntervalMatrix& a, NORMTYPE n );
00397
00398 IntervalMatrix& diagonalScale( IntervalMatrix& d, bool side = true );
00399 friend IntervalMatrix& diagonalScale( IntervalMatrix& a, IntervalMatrix& d,
00400 bool side = true );
00401 static void diagonalScale( IntervalMatrix& a, IntervalMatrix& d,
00402 IntervalMatrix& result, bool side = true );
00403
00404 friend IntervalMatrix& diagonalScale( IntervalMatrix& a, IntervalMatrix& b,
00405 IntervalMatrix& d );
00406 static void diagonalScale( IntervalMatrix& a, IntervalMatrix& b,
00407 IntervalMatrix& d, IntervalMatrix& result );
00408
00409 friend IntervalMatrix& diagonalScale2( IntervalMatrix& a, IntervalMatrix& d1,
00410 IntervalMatrix& d2 );
00411 static void diagonalScale2( IntervalMatrix& a, IntervalMatrix& d1,
00412 IntervalMatrix& d2, IntervalMatrix& result );
00413
00414 friend IntervalMatrix& triangularMult( IntervalMatrix& a, IntervalMatrix& t,
00415 const Interval& alpha = IntervalMatrix::one,
00416 TRIMAT_TYPE tt = LTR_TOP, bool side = false );
00417 static void triangularMult( IntervalMatrix& a, IntervalMatrix& t,
00418 IntervalMatrix& result,
00419 const Interval& alpha = IntervalMatrix::one,
00420 TRIMAT_TYPE tt = LTR_TOP, bool side = false );
00421
00422 friend IntervalMatrix& scaledAccumulation( IntervalMatrix& a, IntervalMatrix& b,
00423 const Interval& alpha = IntervalMatrix::one,
00424 const Interval& beta = IntervalMatrix::one );
00425 static void scaledAccumulation( IntervalMatrix& a, IntervalMatrix& b,
00426 IntervalMatrix& result,
00427 const Interval& alpha = IntervalMatrix::one,
00428 const Interval& beta = IntervalMatrix::one );
00429
00430 friend IntervalMatrix& accTranspose( IntervalMatrix& a, IntervalMatrix& b,
00431 const Interval& alpha = IntervalMatrix::one,
00432 const Interval& beta = IntervalMatrix::one );
00433 static void accTranspose( IntervalMatrix& a, IntervalMatrix& b,
00434 IntervalMatrix& result,
00435 const Interval& alpha = IntervalMatrix::one,
00436 const Interval& beta = IntervalMatrix::one );
00437
00438 friend IntervalVector& scaledVectorMult( IntervalMatrix& a,
00439 IntervalVector& x,
00440 IntervalVector& y,
00441 const Interval& alpha = IntervalMatrix::one,
00442 const Interval& beta = IntervalMatrix::one,
00443 bool transpose_a = false );
00444 static void scaledVectorMult( IntervalMatrix& a,
00445 IntervalVector& x,
00446 IntervalVector& y,
00447 IntervalVector& result,
00448 const Interval& alpha = IntervalMatrix::one,
00449 const Interval& beta = IntervalMatrix::one,
00450 bool transpose_a = false );
00451
00452 friend IntervalMatrix& scaledMatrixMult( IntervalMatrix& a,
00453 IntervalMatrix& b,
00454 IntervalMatrix& c,
00455 bool a_transpose = false,
00456 bool b_transpose = false,
00457 const Interval& alpha = IntervalMatrix::one,
00458 const Interval& beta = IntervalMatrix::one );
00459 static void scaledMatrixMult( IntervalMatrix& a,
00460 IntervalMatrix& b,
00461 IntervalMatrix& c,
00462 IntervalMatrix& result,
00463 bool a_transpose = false,
00464 bool b_transpose = false,
00465 const Interval& alpha = IntervalMatrix::one,
00466 const Interval& beta = IntervalMatrix::one );
00467
00468 friend IntervalMatrix& rankOneUpdate( IntervalMatrix& a,
00469 IntervalVector& x,
00470 IntervalVector& y,
00471 const Interval& alpha = IntervalMatrix::one,
00472 const Interval& beta = IntervalMatrix::one );
00473 static void rankOneUpdate( IntervalMatrix& a,
00474 IntervalVector& x,
00475 IntervalVector& y,
00476 IntervalMatrix& result,
00477 const Interval& alpha = IntervalMatrix::one,
00478 const Interval& beta = IntervalMatrix::one );
00479
00480 IntervalMatrix& permute( int p[] );
00481 friend IntervalMatrix& permute( IntervalMatrix& a, int p[] );
00482 static void permute( IntervalMatrix& a, IntervalMatrix& result, int p[] );
00483
00484
00485 bool encloses( IntervalMatrix& y );
00486 friend bool encloses( IntervalMatrix& v, IntervalMatrix& y );
00487
00488 bool interior( IntervalMatrix& y );
00489 friend bool interior( IntervalMatrix& v, IntervalMatrix& y );
00490
00491 bool disjoint( IntervalMatrix& y );
00492 friend bool disjoint( IntervalMatrix& v, IntervalMatrix& y );
00493
00494
00495 IntervalMatrix& operator&( IntervalMatrix& y );
00496 friend IntervalMatrix& intervalIntersection( IntervalMatrix& v,
00497 IntervalMatrix& y );
00498 static void intervalIntersection( IntervalMatrix& v, IntervalMatrix& y,
00499 IntervalMatrix& result );
00500
00501
00502 IntervalMatrix& operator|( IntervalMatrix& y );
00503 friend IntervalMatrix& intervalUnion( IntervalMatrix& v,
00504 IntervalMatrix& y );
00505 static void intervalUnion( IntervalMatrix& v, IntervalMatrix& y,
00506 IntervalMatrix& result );
00507
00508 IntervalMatrix& intervalHull( IntervalMatrix& y );
00509 friend IntervalMatrix& intervalHull( IntervalMatrix& v, IntervalMatrix& y );
00510 static void intervalHull( IntervalMatrix& v, IntervalMatrix& y,
00511 IntervalMatrix& result );
00512
00513 bool isEmpty();
00514 friend bool isEmpty( IntervalMatrix& v );
00515
00516 IntervalMatrix& lowerBounds();
00517 friend IntervalMatrix& lowerBounds( IntervalMatrix &v );
00518 static void lowerBounds( IntervalMatrix &v, IntervalMatrix& result );
00519
00520 IntervalMatrix& upperBounds();
00521 friend IntervalMatrix& upperBounds( IntervalMatrix &v );
00522 static void upperBounds( IntervalMatrix &v, IntervalMatrix& result );
00523
00524 IntervalMatrix& midpoint();
00525 friend IntervalMatrix& midpoint( IntervalMatrix &v );
00526 static void midpoint( IntervalMatrix &v, IntervalMatrix& result );
00527
00528 IntervalMatrix& widths();
00529 friend IntervalMatrix& widths( IntervalMatrix &v );
00530 static void widths( IntervalMatrix &v, IntervalMatrix& result );
00531
00532
00533 void lu( IntervalMatrix& lu, int* p );
00534 friend void lu( IntervalMatrix& a, IntervalMatrix& lu, int* p );
00535 static void lu( IntervalMatrix& a, IntervalMatrix& lu, int* p );
00536
00537
00538 IntervalVector& solveLU( IntervalVector& b, int p[] );
00539 friend IntervalVector& solveLU( IntervalVector& b, IntervalMatrix& lu, int p[] );
00540 static void solveLU( IntervalVector& b, IntervalMatrix& lu, int p[], IntervalVector& result );
00541
00542 IntervalMatrix& inverse();
00543 friend IntervalMatrix& inverse( IntervalMatrix& a );
00544 static void inverse( IntervalMatrix& a, IntervalMatrix& result );
00545
00546 IntervalVector& toVector( int index, bool row = true );
00547 void fromVector( IntervalVector& v, int index, bool row = true );
00548
00549
00550 static IntervalMatrix& getResult( IntervalMatrix* a, IntervalMatrix* b,
00551 bool force_new = false );
00552 static void freeResult( IntervalMatrix* a, IntervalMatrix* b,
00553 IntervalMatrix& r );
00554 static void freeResult( IntervalMatrix* a, IntervalMatrix* b = NULL );
00555 };
00556
00577 template <unsigned int r, unsigned int c>
00578 class IntervalMatrixT : public IntervalMatrix
00579 {
00580 public:
00581 IntervalMatrixT() : IntervalMatrix( r, c ){ }
00582 IntervalMatrixT( const IntervalMatrix& m ) {
00583 IntervalMatrixT::operator=( m );
00584 }
00585 IntervalMatrixT( const IntervalMatrixT& m ) {
00586 IntervalMatrixT::operator=( m );
00587 }
00588
00599
00613
00632
00795
00812
00839
00864
00870
00959