IntBLAS.h

Go to the documentation of this file.
00001 //IntBLAS: Implementation of the interval BLAS standard.
00002 //Copyright (C) 2005  Micheal Nooner
00003 //
00004 //This library is free software; you can redistribute it and/or
00005 //modify it under the terms of the GNU Lesser General Public
00006 //License as published by the Free Software Foundation; either
00007 //version 2.1 of the License, or (at your option) any later version.
00008 //
00009 //This library is distributed in the hope that it will be useful,
00010 //but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 //MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012 //Lesser General Public License for more details.
00013 //
00014 //You should have received a copy of the GNU Lesser General Public
00015 //License along with this library; if not, write to the Free Software
00016 //Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
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 //Forward declarations.
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   //Error operations
00160   static INTERVAL_EXCEPTION getLastError();
00161   static INTERVAL_EXCEPTION setError( int code, const char* msg = NULL );
00162 
00163   //Gets
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   //Utility operators
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   //Mathematics operators
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; //Power
00203   Interval operator*( double r ) const; //Scalar Multiplication
00204   friend Interval operator*( double r, const Interval& b ) { return (b * r); }
00205   Interval operator+( double r ) const; //Scalar Addition
00206   friend Interval operator+( double r, const Interval& b ) { return (b + r); }
00207   Interval operator-( double r ) const; //Scalar Subtraction
00208   friend Interval operator-( double r, const Interval& b );
00209   Interval operator/( double r ) const; //Scalar Division
00210   friend Interval operator/( double r, const Interval& b );
00211   Interval operator-() const; //Unary negation
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   //Set operations
00230   Interval operator&( const Interval& b ) const; //Intersection
00231   Interval operator|( const Interval& b ) const; //Union
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   //Set tests
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   //Trig functions
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   //Overloads of some cmath functions
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; //The interval [1.0,1.0]
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   //Utility methods
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   //Utility operators
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   //Mathematics operators
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   //Set operations
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   //Intersection
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   //Union
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   //LU Decomposition
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   //Solve equations using LU matrix
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   //Internal Memory Management Methods
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 

Generated on Wed Apr 26 16:12:03 2006 for IntBLAS by  doxygen 1.4.4