00001 #include <stdio.h>
00002 #include <math.h>
00003 #include <string.h>
00004 #include "bmatrix.h"
00005 #include "matrix.h"
00006 #include "vector.h"
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 bmatrix::bmatrix(long m, long n)
00018 {
00019 bmatrix::m = m;
00020 bmatrix::n = n;
00021 a = new matrix[m*n];
00022 memset (a,0,m*n*sizeof(matrix));
00023 if (a == NULL)
00024 print_err("allocating memory for bmatrix", __FILE__, __LINE__, __func__);
00025
00026 #ifdef DEBUG_BMATRIX
00027 Acbm++;
00028 #endif
00029 }
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 bmatrix::bmatrix (const bmatrix &mat)
00040 {
00041 long i, j;
00042 if ((mat.m != m) || (mat.n != n) || (mat.a == NULL))
00043 {
00044 m = mat.m;
00045 n = mat.n;
00046 a = new matrix[m*n];
00047 if (a == NULL)
00048 print_err("allocating memory for bmatrix", __FILE__, __LINE__, __func__);
00049
00050 }
00051 for (i = 0; i < m; i++)
00052 for (j = 0; j < n; j++)
00053 a[i*n+j] = mat.a[i*n+j];
00054 if (mat.row.n)
00055 row = mat.row;
00056 if (mat.col.n)
00057 col = mat.col;
00058 #ifdef DEBUG_BMATRIX
00059 print_warning("copy constructor of bmatrix is called.\n"
00060 "Please check your code and make sure you want use copy constructor.",__FILE__, __LINE__, __func__);
00061 Acbm++;
00062 #endif
00063 }
00064
00065
00066
00067
00068
00069
00070
00071
00072 bmatrix::~bmatrix()
00073 {
00074 m = 0;
00075 n = 0;
00076 #ifdef DEBUG_BMATRIX
00077 if (a != NULL)
00078 Acbm--;
00079 #endif
00080 delete [] a;
00081 a = NULL;
00082 }
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094 long bmatrix::gen_indices()
00095 {
00096 long i, j, tid, d;
00097 allocv(m, row);
00098 allocv(n, col);
00099 for (i=0; i < m; i++)
00100 {
00101 for(j=0; j < n; j++)
00102 {
00103 if (j == 0)
00104 {
00105 row[i] = a[i*n+j].m;
00106 continue;
00107 }
00108 if (row[i] != a[i*n+j].m)
00109 {
00110 print_err("wrong number of rows in matrices in %ld. row of block matrix\n", __FILE__, __LINE__, __func__, i+1);
00111 return 1;
00112 }
00113 }
00114 }
00115 for (i=0; i < m; i++)
00116 {
00117 for(j=0; j < n; j++)
00118 {
00119 if (i == 0)
00120 {
00121 col[j] = a[i*n+j].n;
00122 continue;
00123 }
00124 if (col[j] != a[i*n+j].n)
00125 {
00126 print_err("wrong number of columns in matrices in %ld.column of block matrix", __FILE__, __LINE__, __func__, j+1);
00127 return 2;
00128 }
00129 }
00130 }
00131 tid = 0;
00132 for(i=0; i < m; i++)
00133 {
00134 d = row[i];
00135 row[i] = tid;
00136 tid += d;;
00137 }
00138 tm = tid;
00139 tid = 0;
00140 for(i=0; i < n; i++)
00141 {
00142 d = col[i];
00143 col[i] = tid;
00144 tid += d;;
00145 }
00146 tn = tid;
00147 return 0;
00148 }
00149
00150
00151
00152
00153
00154
00155
00156
00157 long bmatrix::give_totm() const
00158 {
00159 long i, totm = 0;
00160 for (i=0; i < m; i++)
00161 totm += a[i*n].m;
00162 return totm;
00163 }
00164
00165
00166
00167
00168
00169
00170
00171
00172 long bmatrix::give_totn() const
00173 {
00174 long i, totn = 0;
00175 for (i=0; i < n; i++)
00176 totn += a[i].n;
00177 return totn;
00178 }
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193 void bmatrix::give_totid(long i, long j, long k, long l, long &ti, long &tj) const
00194 {
00195 #ifdef DEBUG_BMATRIX
00196 if ((i >= m) || (j >=n) || (a[i][j].m >= k) || (a[i][j]>=l) ||
00197 (i < 0) || (j < 0) || (k < 0) || (l < 0))
00198 {
00199 print_err("required indices are out of range", __FILE__, __LINE__, __func__);
00200 abort();
00201 }
00202 #endif
00203 if (row.n && col.n)
00204 {
00205 ti = row[i]+k;
00206 tj = col[j]+l;
00207 return;
00208 }
00209
00210 long ii;
00211 ti = 0;
00212 for(ii=0; ii<i; ii++)
00213 ti += a[ii*n].m;
00214 ti += k;
00215 tj = 0;
00216 for(ii=0; ii<j; ii++)
00217 tj += a[ii].n;
00218 tj += l;
00219 }
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235 void bmatrix::give_locid(long ti, long tj, long &i, long &j, long &k, long &l) const
00236 {
00237 #ifdef DEBUG_BMATRIX
00238 if (tm == 0)
00239 tm = give_totm();
00240 if (tn == 0)
00241 tn = give_totn();
00242 if ((ti >= tm) || (tj >= tn) || (ti < 0) || (tj < 0))
00243 {
00244 print_err("required indices are out of range", __FILE__, __LINE__, __func__);
00245 abort();
00246 }
00247 #endif
00248 if (row.n && col.n)
00249 {
00250 for(i = 0; i < m; i++)
00251 if (ti < row[i]) break;
00252 for(j = 0; j < n; j++)
00253 if (tj < col[j]) break;
00254 if (i > 0)
00255 k = ti - row[i-1];
00256 else
00257 k = ti;
00258 if (j > 0)
00259 l = tj - col[j-1];
00260 else
00261 l = tj;
00262 return;
00263 }
00264
00265 long d;
00266 d = 0;
00267 for(i=0; i < m; i++)
00268 {
00269 if (ti < d) break;
00270 d += a[i*n].m;
00271 }
00272 if (i > 0)
00273 k = ti - d;
00274 else
00275 k = ti;
00276
00277 d = 0;
00278 for(j=0; j < n; j++)
00279 {
00280 if (tj < d) break;
00281 d += a[j].n;
00282 }
00283 if (j > 0)
00284 l = tj - d;
00285 else
00286 l = tj;
00287 }
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 long allocm(long m, long n, bmatrix &mat)
00301 {
00302 mat.m = m;
00303 mat.n = n;
00304 mat.a = new matrix[m*n];
00305 memset (mat.a,0,m*n*sizeof(matrix));
00306 if (mat.a == NULL)
00307 {
00308 print_err("allocating memory for bmatrix", __FILE__, __LINE__, __func__);
00309 return (1);
00310 }
00311 #ifdef DEBUG_BMATRIX
00312 Acbm++;
00313 #endif
00314 return (0);
00315 }
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 long copym(const bmatrix &src, bmatrix &dest)
00336 {
00337 long i, j;
00338 if ((src.m != dest.m) || (src.n != dest.n))
00339 {
00340 print_err("copying bmatrix - incompatible size of matrices", __FILE__, __LINE__, __func__);
00341 return (1);
00342 }
00343 for (i = 0; i < src.m; i++)
00344 for (j = 0; j < src.n; j++)
00345 dest[i][j] = src[i][j];
00346
00347 if (src.row.n)
00348 dest.row = src.row;
00349 if (src.col.n)
00350 dest.col = src.col;
00351
00352 return (0);
00353 }
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366 long fillm(double c, bmatrix &mat)
00367 {
00368 long i, j, k, l;
00369 for (i = 0; i < mat.m; i++)
00370 for (j = 0; j < mat.n; j++)
00371 for (k = 0; k < mat[i][j].m; k++)
00372 for (l = 0; l < mat[i][j].n; l++)
00373 mat[i][j][k][l] = c;
00374 return (0);
00375 }
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 long destrm(bmatrix &mat)
00390 {
00391 #ifdef DEBUG_BMATRIX
00392 if (mat.a != NULL)
00393 Acbm--;
00394 #endif
00395 mat.m = 0;
00396 mat.n = 0;
00397 delete [] mat.a;
00398 mat.a = NULL;
00399 return (0);
00400 }
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421 long addm(const bmatrix &a, const bmatrix &b, bmatrix &c)
00422 {
00423 long i,j,k,l;
00424 long ti, tj;
00425 long ai, aj, ak, al;
00426 long bi, bj, bk, bl;
00427
00428 if ((a.tm != b.tm) || (a.tn != b.tn) ||
00429 (c.tm != b.tm) || (c.tn != b.tn))
00430 {
00431 print_err("adding bmatrix - incompatible size", __FILE__, __LINE__, __func__);
00432 return (1);
00433 }
00434 for (i = 0; i < c.m; i++)
00435 {
00436 for (j = 0; j < c.n; j++)
00437 {
00438 for (k = 0; k < c[i][j].m; k++)
00439 {
00440 for (l = 0; l < c[i][j].n; l++)
00441 {
00442 c.give_totid(i, j, k, l, ti, tj);
00443 a.give_locid(ti, tj, ai, aj, ak, al);
00444 b.give_locid(ti, tj, bi, bj, bk, bl);
00445 c[i][j][k][l] = a[ai][aj][ak][al] + b[bi][bj][bk][bl];
00446 }
00447 }
00448 }
00449 }
00450 return(0);
00451 }
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471 long subm(const bmatrix &a, const bmatrix &b, bmatrix &c)
00472 {
00473 long i,j,k,l;
00474 long ti, tj;
00475 long ai, aj, ak, al;
00476 long bi, bj, bk, bl;
00477
00478 if ((a.tm != b.tm) || (a.tn != b.tn) ||
00479 (c.tm != b.tm) || (c.tn != b.tn))
00480 {
00481 print_err("adding bmatrix - incompatible size", __FILE__, __LINE__, __func__);
00482 return (1);
00483 }
00484 for (i = 0; i < c.m; i++)
00485 {
00486 for (j = 0; j < c.n; j++)
00487 {
00488 for (k = 0; k < c[i][j].m; k++)
00489 {
00490 for (l = 0; l < c[i][j].n; l++)
00491 {
00492 c.give_totid(i, j, k, l, ti, tj);
00493 a.give_locid(ti, tj, ai, aj, ak, al);
00494 b.give_locid(ti, tj, bi, bj, bk, bl);
00495 c[i][j][k][l] = a[ai][aj][ak][al] + b[bi][bj][bk][bl];
00496 }
00497 }
00498 }
00499 }
00500 return(0);
00501 }
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526 long mxm(const bmatrix &a, const bmatrix &b, bmatrix &c)
00527 {
00528
00529 long i,j,k;
00530 long ai, aj, ak, al;
00531 long bi, bj, bk, bl;
00532 long ci, cj, ck, cl;
00533
00534 if ((a.tn != b.tm) ||
00535 (c.tm != a.tm) || (c.tn != b.tn))
00536 {
00537 print_err("multiplying matrix - incompatible size", __FILE__, __LINE__, __func__);
00538 return (1);
00539 }
00540 for (i=0; i < a.tm; i++)
00541 for (j=0; j < b.tn; j++)
00542 {
00543 c.give_locid(i, j, ci, cj, ck, cl);
00544 c[ci][cj][ck][cl] = 0.0;
00545 for (k=0; k < a.tn; k++)
00546 {
00547 a.give_locid(i, k, ai, aj, ak, al);
00548 b.give_locid(k, j, bi, bj, bk, bl);
00549 c[ci][cj][ck][cl] += a[ai][aj][ak][al] * b[bi][bj][bk][bl];
00550 }
00551 }
00552 return(0);
00553 }
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569 long cmulm(double c, bmatrix &a)
00570 {
00571
00572 long i, j, k, l;
00573 for (i = 0; i < a.m; i++)
00574 for (j = 0; j < a.n; j++)
00575 for (k = 0; k < a[i][j].m; k++)
00576 for (l = 0; l < a[i][j].n; l++)
00577 a[i][j][k][l] *= c;
00578
00579 return (0);
00580 }