00001 #ifndef DEF_MATRIX
00002 #define DEF_MATRIX
00003 #include <assert.h>
00004 #include "vector.h"
00005 #include "gfunct.h"
00006
00007
00008
00009 struct matrix
00010
00011
00012
00013 {
00014 long m;
00015 long n;
00016 double *a;
00017
00018 matrix() {m = n = 0L; a = NULL;};
00019 matrix(long m, long n);
00020 matrix(const matrix &mat);
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 double* operator [] (long i) const
00031 {
00032 #ifdef DEBUG_MATRIX
00033 if (i >= m)
00034 print_err("matrix row index %ld is out of range <0,%ld>", __FILE__, __LINE__, __func__, i, m-1);
00035
00036 assert(i < m);
00037 #endif
00038 return (&a[i * n]);
00039 };
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051 double& operator () (long i, long j) const
00052 {
00053 #ifdef DEBUG_MATRIX
00054 if (i >= m)
00055 print_err("matrix row index %ld is out of range <0,%ld>", __FILE__, __LINE__, __func__, i, m-1);
00056 if (j >= n)
00057 print_err("matrix column index %ld is out of range <0,%ld>", __FILE__, __LINE__, __func__, j, n-1);
00058
00059 assert((i < m) && (j < n));
00060 #endif
00061 return (a[i * n + j]);
00062 };
00063
00064 ~matrix();
00065 };
00066
00067
00068 struct imatrix
00069
00070
00071 {
00072 long m;
00073 long n;
00074 long *a;
00075
00076 imatrix() {m = n = 0L; a = NULL;};
00077 imatrix(long m, long n);
00078 imatrix(const imatrix &mat);
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088 long* operator [] (long i)
00089 {
00090 #ifdef DEBUG_MATRIX
00091 if (i >= m)
00092 print_err("imatrix row index %ld is out of range <0,%ld>", __FILE__, __LINE__, __func__, i, m-1);
00093
00094 assert(i < m);
00095 #endif
00096 return (&a[i * n]);
00097 };
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 long& operator () (long i, long j)
00109 {
00110 #ifdef DEBUG_MATRIX
00111 if (i >= m)
00112 print_err("imatrix row index %ld is out of range <0,%ld>", __FILE__, __LINE__, __func__, i, m-1);
00113 if (j >= n)
00114 print_err("imatrix column index %ld is out of range <0,%ld>", __FILE__, __LINE__, __func__, j, n-1);
00115
00116 assert((i < m) && (j < n));
00117 #endif
00118 return (a[i * n + j]);
00119 };
00120
00121 ~imatrix();
00122 };
00123
00124
00125 struct gfmatrix
00126
00127
00128
00129 {
00130 long m;
00131 long n;
00132 gfunct *a;
00133
00134 gfmatrix() {m = n = 0L; a = NULL;};
00135 gfmatrix(long m, long n);
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 gfunct& operator () (long i, long j) const
00147 {
00148 #ifdef DEBUG_MATRIX
00149 if (i >= m)
00150 print_err("matrix row index %ld is out of range <0,%ld>", __FILE__, __LINE__, __func__, i, m-1);
00151 if (j >= n)
00152 print_err("matrix column index %ld is out of range <0,%ld>", __FILE__, __LINE__, __func__, j, n-1);
00153
00154 assert((i < m) && (j < n));
00155 #endif
00156 return (a[i*n+j]);
00157 };
00158
00159
00160 void evaluate(vector &p, char **namevar, matrix &res);
00161
00162 void evaluate_t(vector &p, char **namevar, matrix &res);
00163
00164 ~gfmatrix();
00165 };
00166
00167
00168
00169 long allocm(long m, long n, matrix &mat);
00170
00171 long allocm(long m, long n, imatrix &mat);
00172
00173 long allocm(long m, long n, double **&mat);
00174
00175 long allocm(long m, long n, gfmatrix &mat);
00176
00177 long allocim(long n, matrix &mat);
00178
00179 long allocim(long n, gfmatrix &mat);
00180
00181
00182 long reallocm(long m, long n, matrix &mat);
00183
00184 long reallocm(long m, long n, imatrix &mat);
00185
00186
00187 long copym(const matrix &src, matrix &dest);
00188
00189
00190 long fillm(double c, matrix &mat);
00191
00192 long fillrow(double c, long i, matrix &mat);
00193
00194 long fillcol(double c, long i, matrix &mat);
00195
00196 long identm(matrix &mat);
00197
00198
00199 long destrm(matrix &mat);
00200
00201 long destrm(double **&mat, long m);
00202
00203 long destrm(long **&mat, long m);
00204
00205
00206 long addm(const matrix &a, const matrix &b, matrix &c);
00207
00208
00209 long subm(const matrix &a, const matrix &b, matrix &c);
00210
00211
00212 void tensprd (vector &a,vector &b,matrix &c);
00213
00214
00215 long mxm(const matrix &a, const matrix &b, matrix &c);
00216
00217 void mxm(const double *a, const double *b, double *c, long l, long m, long n);
00218
00219
00220 long mxmt(const matrix &a, const matrix &b, matrix &c);
00221
00222 void mxmt(const double *a, const double *b, double *c, long l, long m, long n);
00223
00224
00225 long mtxm(const matrix &a, const matrix &b, matrix &c);
00226
00227 void mtxm(const double *a, const double *b, double *c, long l, long m, long n);
00228
00229
00230 void mtxmccr(const double *a, const double *b, double *c, long l, long m, long n);
00231
00232
00233 long cmulm(const double c, const matrix &a, matrix &b);
00234
00235 long cmulm(double c, matrix &a);
00236
00237
00238 long tranm(const matrix &a, matrix &at);
00239
00240 long tranm(matrix &a);
00241
00242
00243 long vxm(const vector &u, const matrix &a, vector &v);
00244
00245 long vxmt(const vector &u, const matrix &a, vector &v);
00246
00247
00248 long mxv(const matrix &a, const vector &u, vector &v);
00249
00250 void mxv(const double *a, const double *b, double *c, long m, long n);
00251
00252
00253 void mxvc(const double *a, const double *b, double *c, long m, long n);
00254
00255 long mixv(const matrix &a, const vector &u, double &vi, long i);
00256
00257
00258 long mtxv(const matrix &a, const vector &u, vector &v);
00259
00260 void mtxv(const double *a, const double *b, double *c, long m, long n);
00261
00262 void mtvc(const double *a, const double *b, double *c, long m, long n);
00263
00264
00265 long vxv(const vector &u, const vector &v, matrix &a);
00266
00267
00268 long vxmxv (const vector &v, const matrix &m, double &answer);
00269
00270 long vxmxv (const double *v, const double *m, long dim, double &answer);
00271
00272
00273
00274 long gause(const matrix &a, matrix &b,double zero);
00275
00276 long gause(const matrix &a, matrix &b, vector &r,double zero);
00277
00278 long gause(const matrix &a, matrix &b, vector &r, vector &sol,double zero);
00279
00280
00281 long invm(const matrix &a, matrix &b,double zero);
00282
00283
00284 long detm(const matrix &a, double &det);
00285
00286
00287 long readm(XFILE *in, matrix &a);
00288
00289 long readm(XFILE *in, gfmatrix &a);
00290
00291
00292 long printm(const matrix &a, FILE *out = stdout, int prec = 3, int width = 11);
00293
00294 long printm(FILE *out, const matrix &a);
00295
00296 long printm(FILE *out, const gfmatrix &a);
00297
00298
00299 long mswapc(matrix &mat, long i, long j);
00300
00301 long mswapr(matrix &mat, long i, long j);
00302
00303
00304 long condense_matrix(matrix &sm, ivector &cu);
00305
00306 long condense_vector(matrix &sm, vector &nf, ivector &cu);
00307
00308
00309 void bdbj (double *k,const double *b,const double *d,double jac,long m,long n);
00310
00311 void bdbjac (matrix &d,const matrix &a,const matrix &b,const matrix &c,double jac);
00312
00313
00314 void nnj (double *m,const double *n,double jac,long mm,long nn);
00315
00316 void nnjac (matrix &d,const matrix &a,const matrix &b,double jac);
00317
00318
00319
00320 void lgvectortransf (vector &g,const vector &l,const matrix &tmat);
00321
00322 void glvectortransf (const vector &g,vector &l,const matrix &tmat);
00323
00324 void glmatrixtransf (matrix &a,const matrix &tmat);
00325
00326 void lgmatrixtransf (matrix &a,const matrix &tmat);
00327
00328
00329 void lgmatrixtransfblock (matrix &a, const matrix &t);
00330
00331 void glvectortransfblock (vector &v, const matrix &t);
00332
00333 void lgvectortransfblock (vector &v, const matrix &t);
00334
00335
00336 void locglob (double *gv, const double *lv,const long *cn,long n);
00337
00338 void globloc (const double *gv, double *lv, const long *cn,long n);
00339
00340
00341 void locvecmat(double *mat, const double *vect, const long *cn, long ci, long m, long n);
00342
00343 void mat_localize (matrix &gm, matrix &lm, long *rcn, long *ccn);
00344
00345
00346 double det2d (const double *x,const double *y);
00347
00348 double det3d (const double *x,const double *y,const double *z);
00349
00350
00351 void princ_val (matrix &v, vector &pval, matrix &pvect, long ni, double err, double zero, long n, long normalize);
00352
00353
00354 void gemp (double *a,double *x,double *y,long n,long m,double limit,long pivot);
00355
00356
00357 void lu_full (double *a,double *x,double *y,long n,double zero,long tc);
00358
00359
00360 void matassem_lsm (double *lsm,vector &natcoord);
00361
00362
00363 void rhsassem_lsm (double *rhs,vector &natcoord,vector &values);
00364
00365
00366 void solve_lsm (double *lsm,double *lhs,double *rhs,double zero,long n,long m);
00367
00368
00369 void extractm (matrix &a,matrix &b,long fi, long ncomp);
00370
00371 void extractblock (double *a,double *b,long n,long nr,long nc,long fri,long fci);
00372
00373 #endif