00001
00002 #ifndef __JBXL_CPP_MATRIX_H_
00003 #define __JBXL_CPP_MATRIX_H_
00004
00005
00006
00015 #include "tools++.h"
00016
00017 #include <math.h>
00018 #include "Vector.h"
00019
00020
00021
00022
00023 namespace jbxl {
00024
00025
00034 template <typename T=double> class Matrix
00035 {
00036 public:
00037 int n;
00038 int r;
00039 int d;
00040 int* sz;
00041 T* mx;
00042 T err;
00043
00044 public:
00045 Matrix() { init();}
00046 Matrix(int nn, ...);
00047 virtual ~Matrix() {}
00048
00049 void init() { n = r = d = 0; err = (T)0; sz = NULL; mx = NULL;}
00050 void init(int nn, ...);
00051
00052 void getm(int nn, int* size);
00053 T& element(int i, ...);
00054 void clear(T v=T(0)) { for(int i=0;i<r;i++) mx[i]=v;}
00055 void dup(Matrix<T> a) { *this = dup_Matrix(a);}
00056
00058 void free() { if(sz!=NULL) ::free(sz); if(mx!=NULL) ::free(mx); init();}
00059 };
00060
00061
00062
00078 template <typename T> Matrix<T>::Matrix(int nn, ...)
00079 {
00080 r = 0;
00081 n = nn;
00082 d = 0;
00083 err = (T)0;
00084 mx = NULL;
00085 sz = (int*)malloc(n*sizeof(int));
00086 if (sz==NULL) return;
00087
00088 va_list argsptr;
00089
00090 va_start(argsptr, nn);
00091 r = 1;
00092 for (int i=0; i<n; i++) {
00093 sz[i] = (int)va_arg(argsptr, int);
00094 r = r*sz[i];
00095 }
00096 va_end(argsptr);
00097
00098 mx = (T*)malloc(r*sizeof(T));
00099 if (mx==NULL) {
00100 ::free(sz);
00101 sz = NULL;
00102 return;
00103 }
00104 for (int i=0; i<r; i++) mx[i] = (T)0;
00105
00106 return;
00107 }
00108
00109
00110
00120 template <typename T> void Matrix<T>::init(int nn, ...)
00121 {
00122 if (sz!=NULL) ::free(sz);
00123 if (mx!=NULL) ::free(mx);
00124
00125 r = 0;
00126 n = nn;
00127 d = 0;
00128 err = (T)0;
00129 mx = NULL;
00130 sz = (int*)malloc(n*sizeof(int));
00131 if (sz==NULL) return;
00132
00133 va_list argsptr;
00134
00135 va_start(argsptr, nn);
00136 r = 1;
00137 for (int i=0; i<n; i++) {
00138 sz[i] = (int)va_arg(argsptr, int);
00139 r = r*sz[i];
00140 }
00141 va_end(argsptr);
00142
00143 mx = (T*)malloc(r*sizeof(T));
00144 if (mx==NULL) {
00145 ::free(sz);
00146 sz = NULL;
00147 return;
00148 }
00149 for (int i=0; i<r; i++) mx[i] = (T)0;
00150
00151 return;
00152 }
00153
00154
00155
00165 template <typename T> void Matrix<T>::getm(int nn, int* size)
00166 {
00167 if (size==NULL) return;
00168 if (sz!=NULL) ::free(sz);
00169 if (mx!=NULL) ::free(mx);
00170
00171 r = 0;
00172 n = nn;
00173 d = 0;
00174 err = (T)0;
00175 mx = NULL;
00176 sz = (int*)malloc(n*sizeof(int));
00177 if (sz==NULL) return;
00178
00179 r = 1;
00180 for (int i=0; i<n; i++) {
00181 sz[i] = size[i];
00182 r = r*sz[i];
00183 }
00184
00185 mx = (T*)malloc(r*sizeof(T));
00186 if (mx==NULL) {
00187 ::free(sz);
00188 sz = NULL;
00189 return;
00190 }
00191 for (int i=0; i<r; i++) mx[i] = (T)0;
00192
00193 return;
00194 }
00195
00196
00197
00213 template <typename T> T& Matrix<T>::element(int i, ...)
00214 {
00215 int* args;
00216 va_list argsptr;
00217
00218 args = (int*)malloc(n*sizeof(int));
00219 if (args==NULL) return err;
00220
00221 va_start(argsptr, i);
00222 args[0] = i;
00223 for (int m=1; m<n; m++) {
00224 args[m] = (int)va_arg(argsptr, int);
00225 }
00226 va_end(argsptr);
00227
00228 int dx = args[0] - 1;
00229 for (int d=1; d<n; d++) dx = dx*sz[d] + args[d] - 1;
00230 ::free(args);
00231
00232 if (dx>=r || dx<0) return err;
00233 return mx[dx];
00234 }
00235
00236
00237
00246 template <typename T> void print_Matrix(FILE* fp, Matrix<T> a)
00247 {
00248 for (int i=0; i<a.r; i++) {
00249 fprintf(fp, " %15lf", (double)a.mx[i]);
00250 if ((i+1)%a.sz[a.n-1]==0) fprintf(fp, "\n");
00251 }
00252 }
00253
00254
00255
00261 template <typename T> Matrix<T> dup_Matrix(Matrix<T> a)
00262 {
00263 Matrix<T> mtx;
00264
00265 if (a.sz!=NULL) {
00266 mtx.getm(a.n, a.sz);
00267 if (a.mx!=NULL && mtx.r==a.r) {
00268 for (int i=0; i<a.r; i++) {
00269 mtx.mx[i] = a.mx[i];
00270 }
00271 }
00272 }
00273
00274 mtx.d = a.d;
00275 mtx.err = (T)(1);
00276
00277 return mtx;
00278 }
00279
00280
00281
00282
00284
00285
00297 template <typename T> Matrix<T> operator * (const Matrix<T> a, const Matrix<T> b)
00298 {
00299 int i, j, k, n, ii, aa, bb;
00300 int *sz, *sa, *sb, *sc, *cx;
00301 T st;
00302 Matrix<T> c;
00303
00304 if (a.mx==NULL || b.mx==NULL) return c;
00305 if (a.sz[a.n-1]!=b.sz[0]) return c;
00306
00307 n = a.n + b.n - 2;
00308 sz = (int*)malloc(n*sizeof(int));
00309 if (sz==NULL) return c;
00310 sa = (int*)malloc(a.n*sizeof(int));
00311 if (sa==NULL) {free(sz); return c;}
00312 sb = (int*)malloc(b.n*sizeof(int));
00313 if (sb==NULL) {free(sz); free(sa); return c;}
00314 sc = (int*)malloc(n*sizeof(int));
00315 if (sc==NULL) {free(sz); free(sa); free(sb); return c;}
00316 cx = (int*)malloc(n*sizeof(int));
00317 if (cx==NULL) {free(sz); free(sa); free(sb); free(sc); return c;}
00318
00319 for (i=0; i<a.n-1; i++) sz[i] = a.sz[i];
00320 for (i=1; i<b.n; i++) sz[a.n-2+i] = b.sz[i];
00321
00322 sa[a.n-1] = sb[b.n-1] = sc[n-1] = 1;
00323 for (i=a.n-2; i>=0; i--) sa[i] = sa[i+1]*a.sz[i+1];
00324 for (i=b.n-2; i>=0; i--) sb[i] = sb[i+1]*b.sz[i+1];
00325 for (i=n-2; i>=0; i--) sc[i] = sc[i+1]*sz[i+1];
00326
00327 c.getm(n, sz);
00328
00329 if (c.sz!=NULL) {
00330 for (i=0; i<c.r; i++) {
00331 ii = i;
00332 for (j=0; j<c.n; j++) {
00333 cx[j] = ii / sc[j];
00334 ii = ii % sc[j];
00335 }
00336 aa = bb = 0;
00337 for (j=0; j<a.n-1; j++) aa = aa + sa[j]*cx[j];
00338 for (j=1; j<b.n; j++) bb = bb + sb[j]*cx[j+a.n-2];
00339
00340 st = (T)0;
00341 for (k=0; k<b.sz[0]; k++) st = st + a.mx[k+aa]*b.mx[bb+sb[0]*k];
00342 c.mx[i] = st;
00343 }
00344 }
00345
00346 free(sz);
00347 free(sa);
00348 free(sb);
00349 free(sc);
00350 free(cx);
00351
00352 return c;
00353 }
00354
00355
00356
00368 template <typename T> Vector<T> operator * (const Matrix<T> a, const Vector<T> v)
00369 {
00370 Vector<T> vct(0.0, 0.0, 0.0, 0.0, -1.0);
00371
00372 Matrix<T> b(1, 3);
00373 b.mx[0] = v.x;
00374 b.mx[1] = v.y;
00375 b.mx[2] = v.z;
00376
00377 Matrix<T> c = a*b;
00378 b.free();
00379
00380 if (c.mx==NULL) return vct;
00381 vct.x = c.mx[0];
00382 vct.y = c.mx[1];
00383 vct.z = c.mx[2];
00384 vct.c = v.c;
00385 c.free();
00386
00387 return vct;
00388 }
00389
00390
00391
00392 template <typename T> inline Matrix<T> operator - (const Matrix<T> a)
00393 {
00394 Matrix<T> c(a.n, a.sz);
00395 if (c.mx!=NULL) for (int i=0; i<a.r; i++) c.mx[i] = -a.mx[i];
00396 return c;
00397 }
00398
00399
00400
00401 template <typename T> inline Matrix<T> operator + (const Matrix<T> a, const Matrix<T> b)
00402 {
00403 Matrix<T> c;
00404 if (!isSameDimension(a, b)) return c;
00405
00406 c.getm(a.n, a.sz);
00407 if (c.mx!=NULL) for (int i=0; i<a.r; i++) c.mx[i] = a.mx[i] + b.mx[i];
00408 return c;
00409 }
00410
00411
00412
00413 template <typename T> inline Matrix<T> operator - (const Matrix<T> a, const Matrix<T> b)
00414 {
00415 Matrix<T> c;
00416 if (!isSameDimension(a, b)) return c;
00417
00418 c.getm(a.n, a.sz);
00419 if (c.mx!=NULL) for (int i=0; i<a.r; i++) c.mx[i] = a.mx[i] - b.mx[i];
00420 return c;
00421 }
00422
00423
00424
00425
00426 template <typename T, typename R> inline Matrix<T> operator * (const R d, const Matrix<T> a)
00427 {
00428 Matrix<T> c(a.n, a.sz);
00429 if (c.mx!=NULL) for (int i=0; i<a.r; i++) c.mx[i] = (T)(d)*a.mx[i];
00430 return c;
00431 }
00432
00433
00434
00435
00436 template <typename T, typename R> inline Matrix<T> operator * (const Matrix<T> a, const R d)
00437 {
00438 Matrix<T> c(a.n, a.sz);
00439 if (c.mx!=NULL) for (int i=0; i<a.r; i++) c.mx[i] = a.mx[i]*(T)d;
00440 return c;
00441 }
00442
00443
00444
00445
00446 template <typename T, typename R> inline Matrix<T> operator / (const Matrix<T> a, const R d)
00447 {
00448 Matrix<T> c(a.n, a.sz);
00449 if (c.mx!=NULL) for (int i=0; i<a.r; i++) c.mx[i] = a.mx[i]/(T)d;
00450 return c;
00451 }
00452
00453
00454
00455 template <typename T> inline bool operator == (const Matrix<T> v1, const Matrix<T> v2)
00456 {
00457 if (!isSameDimension(v1, v2)) return false;
00458 for (int i=0; i<v1.r; i++) if (v1.mx[i]!=v2.mx[i]) return false;
00459 return true;
00460 }
00461
00462
00463
00464 template <typename T> inline bool isSameDimension(const Matrix<T> v1, const Matrix<T> v2)
00465 {
00466 if (v1.n!=v2.n || v1.r!=v2.r) return false;
00467 for (int i=0; i<v1.n; i++) if (v1.sz[i]!=v2.sz[i]) return false;
00468 return true;
00469 }
00470
00471
00472 }
00473
00474
00475 #endif
