00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
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
00026
00027
00028
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
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
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
00068
00069
00070
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
00088
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
00114
00115
00116
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
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
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
00145
00146
00147
00148
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
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
00200
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
00258
00259
00260
00261
00262
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
00274
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 }