00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef SCPPNT_MATOP_H
00022 #define SCPPNT_MATOP_H
00023
00024 #ifdef SCPPNT_NO_DIR_PREFIX
00025 #include "scppnt.h"
00026 #include "slice_pointer.h"
00027 #include "vec.h"
00028 #include "scppnt_error.h"
00029 #else
00030 #include "scppnt/scppnt.h"
00031 #include "scppnt/slice_pointer.h"
00032 #include "scppnt/vec.h"
00033 #include "scppnt/scppnt_error.h"
00034 #endif
00035
00036 namespace SCPPNT
00037 {
00038
00039
00040 template<class MAT> MAT transpose(const MAT &A)
00041 {
00042 Subscript M = A.num_rows();
00043 Subscript N = A.num_columns();
00044
00045 MAT S(N, M);
00046 Subscript i, j;
00047
00048 typename MAT::const_rows_iterator rowsA = A.begin_rows();
00049 typename MAT::columns_iterator columnsS = S.begin_columns();
00050 for (i=M; i--; ++rowsA, ++columnsS)
00051 {
00052 typename MAT::const_row_iterator rowA = *rowsA;
00053 typename MAT::column_iterator columnS = *columnsS;
00054 for (j=N; j--; ++rowA, ++columnS)
00055 *columnS = *rowA;
00056 }
00057 return S;
00058 }
00059
00060
00061 template<class MA, class MB> void matcopy(MA &A, const MB &B)
00062 {
00063 if (A.num_rows() != B.num_rows() || A.num_columns() != B.num_columns())
00064 throw BadDimension("SCPPNT::matcopy");
00065
00066 typename MB::const_rows_iterator brows = B.begin_rows();
00067 typename MA::rows_iterator arows = A.begin_rows();
00068 for (int i = A.num_rows(); i--; ++brows, ++arows)
00069 {
00070 typename MB::const_row_iterator brow = *brows;
00071 typename MA::row_iterator arow = *arows;
00072 for (int j = A.num_columns(); j--; ++arow, ++brow)
00073 {
00074 *arow = *brow;
00075 }
00076 }
00077
00078 }
00079
00080
00081 template<class MA, class MB, class MC> void matadd(const MA &A, const MB &B, MC &C)
00082 {
00083 Subscript nr = A.num_rows();
00084 Subscript nc = A.num_columns();
00085 if (nr != B.num_rows() || nc != B.num_columns() || nr != C.num_rows() || nc != C.num_columns())
00086 throw BadDimension("SCPPNT::matadd(A, B, C)");
00087
00088 typename MB::const_rows_iterator brows = B.begin_rows();
00089 typename MA::const_rows_iterator arows = A.begin_rows();
00090 typename MC::rows_iterator crows = C.begin_rows();
00091 for (int i = nr; i--; ++brows, ++arows, ++crows)
00092 {
00093 typename MB::const_row_iterator brow = *brows;
00094 typename MA::const_row_iterator arow = *arows;
00095 typename MC::row_iterator crow = *crows;
00096 for (int j = nc; j--; ++arow, ++brow, ++crow)
00097 {
00098 *crow = *arow + *brow;
00099 }
00100 }
00101 }
00102
00103
00104 template<class MA, class MB> void matadd(MA &A, const MB &B)
00105 {
00106
00107 if (A.num_rows() != B.num_rows() || A.num_columns() != B.num_columns())
00108 throw BadDimension("SCPPNT::matadd");
00109
00110 typename MB::const_rows_iterator brows = B.begin_rows();
00111 typename MA::rows_iterator arows = A.begin_rows();
00112 for (int i = A.num_rows(); i--; ++brows, ++arows)
00113 {
00114 typename MB::const_row_iterator brow = *brows;
00115 typename MA::row_iterator arow = *arows;
00116 for (int j = A.num_columns(); j--; ++arow, ++brow)
00117 {
00118 *arow += *brow;
00119 }
00120 }
00121
00122 }
00123
00124
00125 template<class MA, class MB, class MC> void matsub(const MA &A, const MB &B, MC &C)
00126 {
00127 Subscript nr = A.num_rows();
00128 Subscript nc = A.num_columns();
00129 if (nr != B.num_rows() || nc != B.num_columns() || nr != C.num_rows() || nc != C.num_columns())
00130 throw BadDimension("SCPPNT::matsub(A, B, C)");
00131
00132 typename MB::const_rows_iterator brows = B.begin_rows();
00133 typename MA::const_rows_iterator arows = A.begin_rows();
00134 typename MC::rows_iterator crows = C.begin_rows();
00135 for (int i = nr; i--; ++brows, ++arows, ++crows)
00136 {
00137 typename MB::const_row_iterator brow = *brows;
00138 typename MA::const_row_iterator arow = *arows;
00139 typename MC::row_iterator crow = *crows;
00140 for (int j = nc; j--; ++arow, ++brow, ++crow)
00141 {
00142 *crow = *arow - *brow;
00143 }
00144 }
00145 }
00146
00147
00148 template<class MA, class MB> void matsub(MA &A, const MB &B)
00149 {
00150
00151 if (A.num_rows() != B.num_rows() || A.num_columns() != B.num_columns())
00152 throw BadDimension("SCPPNT::matsub");
00153
00154 typename MB::const_rows_iterator brows = B.begin_rows();
00155 typename MA::rows_iterator arows = A.begin_rows();
00156 for (int i = A.num_rows(); i--; ++brows, ++arows)
00157 {
00158 typename MB::const_row_iterator brow = *brows;
00159 typename MA::row_iterator arow = *arows;
00160 for (int j = A.num_columns(); j--; ++arow, ++brow)
00161 {
00162 *arow -= *brow;
00163 }
00164 }
00165 }
00166
00167
00168
00169
00170
00171 template<class MR, class MA, class MB> MR mult_element(const MA &A, const MB &B)
00172 {
00173 Subscript M = A.num_rows();
00174 Subscript N = A.num_columns();
00175
00176 if (M != B.num_rows() || N != B.num_columns())
00177 throw BadDimension("SCPPNT::mult_element");
00178
00179 MR tmp(M, N);
00180 Subscript i, j;
00181
00182 typename MA::const_rows_iterator rowsA = A.begin_rows();
00183 typename MB::const_rows_iterator rowsB = B.begin_rows();
00184 typename MR::rows_iterator rowstmp = tmp.begin_rows();
00185 for (i=M; i--; ++rowsA, ++rowsB, ++rowstmp)
00186 {
00187 typename MA::const_row_iterator rowA = *rowsA;
00188 typename MB::const_row_iterator rowB = *rowsB;
00189 typename MR::row_iterator rowtmp = *rowstmp;
00190 for (j=N; j--; ++rowA, ++rowB, ++rowtmp)
00191 *rowtmp = *rowA * *rowB;
00192 }
00193 return tmp;
00194 }
00195
00196
00197 template<class MR, class MA, class MB> MR matmult(const MA &A, const MB &B)
00198 {
00199 Subscript M = A.num_rows();
00200 Subscript N = A.num_columns();
00201 Subscript K = B.num_columns();
00202
00203 if (N != B.num_rows())
00204 throw BadDimension("SCPPNT::matmult(A, B)");
00205
00206 MR tmp(M, K);
00207 typename MR::element_type sum;
00208
00209 typename MR::rows_iterator rowsTmp = tmp.begin_rows();
00210 typename MA::const_rows_iterator rowsA = A.begin_rows();
00211 for (Subscript i=M; i--; ++rowsTmp, ++rowsA)
00212 {
00213 typename MR::row_iterator rowTmp = *rowsTmp;
00214 typename MB::const_columns_iterator columnsB = B.begin_columns();
00215 for (Subscript k=K; k--; ++rowTmp, ++columnsB)
00216 {
00217 typename MA::const_row_iterator rowA = *rowsA;
00218 typename MB::const_column_iterator columnB = *columnsB;
00219 sum = 0;
00220 for (Subscript j=N; j--; ++rowA, ++columnB)
00221 sum += *rowA * *columnB;
00222
00223 *rowTmp = sum;
00224 }
00225 }
00226 return tmp;
00227 }
00228
00229
00230 template<class MC, class MA, class MB> void matmult(MC &C, const MA &A, const MB &B)
00231 {
00232 Subscript M = A.num_rows();
00233 Subscript N = A.num_columns();
00234 Subscript K = B.num_columns();
00235
00236 if (N != B.num_rows())
00237 throw BadDimension("SCPPNT::matmult(C, A, B)");
00238
00239 C.newsize(M, K);
00240
00241 typename MC::element_type sum;
00242
00243 typename MC::rows_iterator rowsC = C.begin_rows();
00244 typename MA::const_rows_iterator rowsA = A.begin_rows();
00245 for (Subscript i=M; i--; ++rowsC, ++rowsA)
00246 {
00247 typename MC::row_iterator rowC = *rowsC;
00248 typename MB::const_columns_iterator columnsB = B.begin_columns();
00249 for (Subscript k=K; k--; ++rowC, ++columnsB)
00250 {
00251 typename MA::const_row_iterator rowA = *rowsA;
00252 typename MB::const_column_iterator columnB = *columnsB;
00253 sum = 0;
00254 for (Subscript j=N; j--; ++rowA, ++columnB)
00255 sum += *rowA * *columnB;
00256
00257 *rowC = sum;
00258 }
00259 }
00260 }
00261
00262
00263 template<class MA, class MB> void matmult_assign(MA &A, const MB &B)
00264 {
00265 Subscript i, j, k;
00266 Subscript M = A.num_rows();
00267 Subscript N = A.num_columns();
00268
00269
00270 if (N != B.num_rows() || B.num_columns() != B.num_rows())
00271 throw BadDimension("SCPPNT::matmult_assign");
00272
00273
00274
00275
00276 typedef typename MA::element_type A_element_type;
00277 A_element_type *temp_vec = new A_element_type[N];
00278
00279 typename MA::rows_iterator irows = A.begin_rows();
00280 for (i = M; i--; ++irows)
00281 {
00282 typename MB::const_columns_iterator icolumns = B.begin_columns();
00283
00284
00285
00286 A_element_type *iv = temp_vec;
00287 typename MB::const_column_iterator icolumn = *icolumns;
00288 typename MA::row_iterator irow = *irows;
00289 A_element_type sum = 0.0;
00290 for (j = N; j--; ++irow, ++icolumn)
00291 {
00292 *iv = *irow;
00293 ++iv;
00294 sum += *irow * *icolumn;
00295 }
00296
00297
00298
00299 irow = *irows;
00300 *irow = sum;
00301
00302
00303 ++irow;
00304 ++icolumns;
00305 for (j=N-1; j--; ++icolumns, ++irow)
00306 {
00307 iv = temp_vec;
00308 icolumn = *icolumns;
00309 sum = 0.0;
00310 for (k=N; k--; ++iv, ++icolumn)
00311 sum += *iv * *icolumn;
00312 *irow = sum;
00313 }
00314 }
00315 delete [] temp_vec;
00316 }
00317
00318
00319
00320
00321
00322
00323
00324
00325 template<class MAT, class IT1, class IT2> void MatTimesVec(const MAT &A, IT1 begin, IT2 result)
00326 {
00327
00328 typedef typename MAT::value_type value_type;
00329
00330 Subscript N = A.num_columns();
00331 Subscript M = A.num_rows();
00332
00333 typename MAT::const_rows_iterator irows = A.begin_rows();
00334 for (Subscript i=M; i--; ++irows, ++result)
00335 {
00336 value_type sum = 0;
00337 typename MAT::const_row_iterator rowi = *irows;
00338 IT1 iv = begin;
00339 for (Subscript j=N; j--; ++iv, ++rowi)
00340 {
00341 sum += *rowi * *iv;
00342 }
00343
00344 *result = sum;
00345 }
00346 }
00347
00348
00349 template<class M, class V> inline V matrix_times_vector(const M &A, const V &x)
00350 {
00351 if (A.num_columns() != x.size())
00352 throw BadDimension("SCPPNT::matrix_times_vector()");
00353
00354 V tmp(A.num_rows());
00355 MatTimesVec(A, x.begin(), tmp.begin());
00356
00357 return tmp;
00358 }
00359
00360 }
00361
00362 #endif
00363