Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

matrix_sub.cpp

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                      Linear Algebra Package
00006  *
00007  *              Basic linear algebra operations, levels 1 & 2
00008  *          Operations on a single row, column, or the diagonal
00009  *                              of a matrix
00010  *
00011  * $Id: matrix_sub.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00012  *
00013  ************************************************************************
00014  */
00015 
00016 #include "LAStreams.h"
00017 #include <math.h>
00018 #include <iostream.h>
00019 
00020 namespace linalg 
00021 {
00022 
00023 /*
00024  *------------------------------------------------------------------------
00025  *              Messing with a single row/column/diag of a matrix
00026  */
00027 
00028                         // Constructing the MatrixColumn object
00029 ConstMatrixColumn::ConstMatrixColumn(const Matrix& m, const int col)
00030         : Matrix::ConstReference(m),
00031           DimSpec(m.q_row_lwb(),m.q_row_upb(),col,col),
00032           col_ptr(m.elements+(col-m.q_col_lwb())*m.q_nrows())
00033 {
00034   if( col > m.q_col_upb() || col < m.q_col_lwb() )
00035     m.info(),
00036     _error("Column #%d is not within the above matrix",col);
00037 }
00038 
00039                         // Constructing the MatrixRow object
00040 ConstMatrixRow::ConstMatrixRow(const Matrix& m, const int row)
00041         : Matrix::ConstReference(m),
00042           DimSpec(row,row,m.q_col_lwb(),m.q_col_upb()),
00043           row_ptr(m.elements+(row-m.q_row_lwb())),
00044           stride(m.q_nrows()),
00045           end_ptr(m.elements+m.nelems)
00046 {
00047   if( row > m.q_row_upb() || row < m.q_row_lwb() )
00048     m.info(),
00049     _error("Row #%d is not within the above matrix",row);
00050   assert( stride > 0 );
00051 }
00052 
00053                         // Constructing the MatrixDiag object
00054 ConstMatrixDiag::ConstMatrixDiag(const Matrix& m)
00055         : Matrix::ConstReference(m),
00056           DimSpec(1,min(m.q_nrows(),m.q_ncols()),
00057                   1,min(m.q_nrows(),m.q_ncols())),
00058           start_ptr(m.elements),
00059           stride(m.q_nrows()+1),
00060           end_ptr(m.elements+m.nelems)
00061 {
00062   assert( stride > 1 );
00063 }
00064 
00065 /*
00066  *------------------------------------------------------------------------
00067  *      Collection-scalar arithmetics: walking with a stride
00068  */
00069 
00070                                 // For every element, do `elem OP value`
00071 #define COMPUTED_VAL_ASSIGNMENT(OP,VALTYPE)                             \
00072                                                                         \
00073 void ElementWiseStride::operator OP (const VALTYPE val)                 \
00074 {                                                                       \
00075   for(register REAL * ep = start_ptr; ep < end_ptr; ep += stride)       \
00076     *ep OP val;                                                         \
00077 }                                                                       \
00078 
00079 COMPUTED_VAL_ASSIGNMENT(=,REAL)
00080 COMPUTED_VAL_ASSIGNMENT(+=,double)
00081 COMPUTED_VAL_ASSIGNMENT(-=,double)
00082 COMPUTED_VAL_ASSIGNMENT(*=,double)
00083 
00084 #undef COMPUTED_VAL_ASSIGNMENT
00085 
00086 
00087                                 // is "element OP val" true for all
00088                                 // elements in the collection?
00089 
00090 #define COMPARISON_WITH_SCALAR(OP)                                      \
00091                                                                         \
00092 bool ElementWiseStrideConst::operator OP (const REAL val) const         \
00093 {                                                                       \
00094   for(register const REAL * ep = start_ptr; ep < end_ptr; ep += stride) \
00095     if( !(*ep OP val) )                                                 \
00096       return false;                                                     \
00097                                                                         \
00098   return true;                                                          \
00099 }                                                                       \
00100 
00101 
00102 COMPARISON_WITH_SCALAR(==)
00103 COMPARISON_WITH_SCALAR(!=)
00104 COMPARISON_WITH_SCALAR(<)
00105 COMPARISON_WITH_SCALAR(<=)
00106 COMPARISON_WITH_SCALAR(>)
00107 COMPARISON_WITH_SCALAR(>=)
00108 
00109 #undef COMPARISON_WITH_SCALAR
00110 
00111 /*
00112  *------------------------------------------------------------------------
00113  *      Apply algebraic functions to all elements of a collection
00114  */
00115 
00116                                 // Take an absolute value of a matrix
00117 void ElementWiseStride::abs(void)
00118 {
00119   for(register REAL * ep = start_ptr; ep < end_ptr; ep += stride)
00120     *ep = fabs(*ep);
00121 }
00122 
00123                                 // Square each element
00124 void ElementWiseStride::sqr(void)
00125 {
00126   for(register REAL * ep = start_ptr; ep < end_ptr; ep += stride )
00127     *ep = *ep * *ep;
00128 }
00129 
00130                                 // Take the square root of all the elements
00131 void ElementWiseStride::sqrt(void)
00132 {
00133   for(register REAL * ep = start_ptr; ep < end_ptr; ep += stride )
00134     if( *ep >= 0 )
00135       *ep = ::sqrt(*ep);
00136     else
00137       _error("%d-th element, %g, is negative. Can't take the square root",
00138              (ep-start_ptr), *ep );
00139 }
00140 
00141 
00142 /*
00143  *------------------------------------------------------------------------
00144  *              Element-wise operations on two groups of elements
00145  */
00146 
00147 
00148                                 // For every element, do `elem OP another.elem`
00149 #define TWO_GROUP_COMP(OP)                                      \
00150                                                                         \
00151 bool ElementWiseStrideConst::operator OP (const ElementWiseStrideConst& another) const \
00152 {                                                                       \
00153   register const REAL * sp = another.start_ptr;                         \
00154   register const REAL * tp = start_ptr;                                 \
00155   for(; tp < end_ptr && sp < another.end_ptr;                           \
00156       tp += stride, sp += another.stride)                               \
00157     if( !(*tp OP *sp) )                                                 \
00158       return false;                                                     \
00159                                                                         \
00160   assure( tp >= end_ptr && sp >= another.end_ptr,                       \
00161     "stride collections have different number of elements" );           \
00162   return true;                                                          \
00163 }                                                                       \
00164 
00165 TWO_GROUP_COMP(==)
00166 TWO_GROUP_COMP(!=)
00167 TWO_GROUP_COMP(<)
00168 TWO_GROUP_COMP(<=)
00169 TWO_GROUP_COMP(>)
00170 TWO_GROUP_COMP(>=)
00171 
00172 #undef TWO_GROUP_COMP
00173 
00174                                 // For every element, do `elem OP another.elem`
00175 #define TWO_GROUP_OP(OP)                                        \
00176                                                                         \
00177 void ElementWiseStride::operator OP (const ElementWiseStrideConst& another) \
00178 {                                                                       \
00179   register const REAL * sp = another.start_ptr;                         \
00180   register REAL * tp = start_ptr;                                       \
00181   for(; tp < end_ptr && sp < another.end_ptr;                           \
00182       tp += stride, sp += another.stride)                               \
00183     *tp OP *sp;                                                         \
00184   assure( tp >= end_ptr && sp >= another.end_ptr,                       \
00185     "stride collections have different number of elements" );           \
00186 }                                                                       \
00187 
00188 TWO_GROUP_OP(=)
00189 TWO_GROUP_OP(+=)
00190 TWO_GROUP_OP(-=)
00191 TWO_GROUP_OP(*=)
00192 TWO_GROUP_OP(/=)
00193 
00194 #undef TWO_GROUP_OP
00195 
00196 
00197 /*
00198  *------------------------------------------------------------------------
00199  *      Reduce a collection or a difference between two collections
00200  *              to a single number: a "norm"
00201  */
00202 
00203 #define REDUCE_SUM(X,VAL) X += (VAL)
00204 #define REDUCE_SUMSQ(X,VAL) X += sqr(VAL)
00205 #define REDUCE_SUMABS(X,VAL) X += fabs(VAL)
00206 #define REDUCE_MAXABS(X,VAL) X = max((REAL)X,fabs(VAL))
00207 
00208 #define REDUCE_ONE(NAME,OP)                                                     \
00209                                                                         \
00210 double ElementWiseStrideConst::NAME (void) const                        \
00211 {                                                                       \
00212   register double norm = 0;                                             \
00213   for(register const REAL * ep = start_ptr; ep < end_ptr; ep += stride) \
00214     OP(norm,*ep);                                                       \
00215   return norm;                                                          \
00216 }                                                                       \
00217 
00218 
00219 
00220 REDUCE_ONE(sum,REDUCE_SUM)
00221 REDUCE_ONE(sum_squares,REDUCE_SUMSQ)
00222 REDUCE_ONE(sum_abs,REDUCE_SUMABS)
00223 REDUCE_ONE(max_abs,REDUCE_MAXABS)
00224 
00225 #undef REDUCE_ONE
00226 
00227 #define REDUCE_DIFF_OF_TWO(NAME,OP)                                     \
00228                                                                         \
00229 double ElementWiseStrideConst::NAME (const ElementWiseStrideConst& another) const       \
00230 {                                                                       \
00231   register double norm = 0;                                             \
00232   register const REAL * sp = another.start_ptr;                         \
00233   register const REAL * tp = start_ptr;                                 \
00234   for(; tp < end_ptr && sp < another.end_ptr;                           \
00235       tp += stride, sp += another.stride)                               \
00236     OP(norm,*tp - *sp);                                                 \
00237                                                                         \
00238   assure( tp >= end_ptr && sp >= another.end_ptr,                       \
00239     "stride collections have different number of elements" );           \
00240   return norm;                                                          \
00241 }                                                                       \
00242 
00243 REDUCE_DIFF_OF_TWO(sum_squares,REDUCE_SUMSQ)
00244 REDUCE_DIFF_OF_TWO(sum_abs,REDUCE_SUMABS)
00245 REDUCE_DIFF_OF_TWO(max_abs,REDUCE_MAXABS)
00246 
00247 #undef REDUCE_DIFF_OF_TWO
00248 
00249 #undef REDUCE_SUM
00250 #undef REDUCE_SUMSQ
00251 #undef REDUCE_SUMABS
00252 #undef REDUCE_MAXABS
00253 
00254 
00255 /*
00256  *------------------------------------------------------------------------
00257  *                 Multiplications with the diagonal matrix
00258  */
00259 
00260                                 // Multiply a matrix by the diagonal
00261                                 // of another matrix
00262                                 // matrix(i,j) *= diag(j)
00263 Matrix& operator *= (Matrix& m, const ConstMatrixDiag& diag)
00264 {
00265   LAStreamOut m_str(m);
00266   LAStrideStreamIn diag_str(diag);
00267 
00268 
00269   if( m.q_ncols() != diag.q_nrows() )
00270     m.info(), cerr << "\n and the diagonal " << diag << endl,
00271    _error("cannot be multiplied");
00272 
00273                         // Each column of m gets multiplied by the corresponding
00274                         // diag(i). Note that m is traversed column-wise
00275   while( !diag_str.eof() )
00276   {
00277     const REAL diag_el = diag_str.get();
00278     for(register int i=0; i<m.q_nrows(); i++)
00279       m_str.get() *= diag_el;
00280   }
00281   assert( m_str.eof() );
00282   return m;
00283 }
00284 }

Generated on Wed Dec 15 21:20:29 2004 for vuVolume by  doxygen 1.3.9.1