00001
00008 #include "matrix.h"
00009
00010
00011
00020 vector unit_vector(vector a)
00021 {
00022 vector c;
00023 double r;
00024
00025 r = a.x*a.x+a.y*a.y+a.z*a.z;
00026 if (Xabs(r)<EPS*EPS) {
00027 c.x = 0.0;
00028 c.y = 0.0;
00029 c.z = 0.0;
00030 c.n = 1.0;
00031 }
00032 else {
00033 r = sqrt(r);
00034 c.x = a.x/r;
00035 c.y = a.y/r;
00036 c.z = a.z/r;
00037 c.n = 1.0;
00038 }
00039 return c;
00040 }
00041
00042
00043
00052 vector unit_ivector(ivector a)
00053 {
00054 vector c;
00055 double r;
00056
00057 r = a.x*a.x+a.y*a.y+a.z*a.z;
00058 if (Xabs(r)<EPS*EPS) {
00059 c.x = 0.0;
00060 c.y = 0.0;
00061 c.z = 0.0;
00062 c.n = 1.0;
00063 }
00064 else {
00065 r = sqrt(r);
00066 c.x = a.x/r;
00067 c.y = a.y/r;
00068 c.z = a.z/r;
00069 c.n = 1.0;
00070 }
00071 return c;
00072 }
00073
00074
00075
00086 vector set_vector(double x, double y, double z)
00087 {
00088 vector c;
00089
00090 c.x = x;
00091 c.y = y;
00092 c.z = z;
00093 c.n = sqrt(x*x + y*y + z*z);
00094 return c;
00095 }
00096
00097
00098
00109 ivector set_ivector(int x, int y, int z)
00110 {
00111 ivector c;
00112
00113 c.x = x;
00114 c.y = y;
00115 c.z = z;
00116 c.n = sqrt((double)x*x + y*y + z*z);
00117 return c;
00118 }
00119
00120
00121
00131 ivector f2ivector(vector a)
00132 {
00133 ivector c;
00134
00135 c.x = (int)(a.x + 0.5);
00136 c.y = (int)(a.y + 0.5);
00137 c.z = (int)(a.z + 0.5);
00138 c.n = sqrt((double)c.x*c.x + c.y*c.y + c.z*c.z);
00139 return c;
00140 }
00141
00142
00143
00152 vector i2vector(ivector a)
00153 {
00154 vector c;
00155
00156 c.x = (double)a.x;
00157 c.y = (double)a.y;
00158 c.z = (double)a.z;
00159 c.n = sqrt(c.x*c.x + c.y*c.y + c.z*c.z);
00160 return c;
00161 }
00162
00163
00164
00175 vector ex_vector(vector a, vector b)
00176 {
00177 vector c;
00178
00179 c.x = a.y*b.z - a.z*b.y;
00180 c.y = a.z*b.x - a.x*b.z;
00181 c.z = a.x*b.y - a.y*b.x;
00182 c.n = sqrt(c.x*c.x+c.y*c.y+c.z*c.z);
00183 return c;
00184 }
00185
00186
00187
00197 matrix make_matrix1(int n)
00198 {
00199 int i;
00200 matrix a;
00201
00202 a.n = a.r = 0;
00203 a.mx = NULL;
00204 a.sz = (int*)malloc(sizeof(int));
00205 if (a.sz==NULL) return a;
00206 a.sz[0] = n;
00207
00208 a.mx = (double*)malloc(n*sizeof(double));
00209 if (a.mx==NULL) {
00210 free(a.sz);
00211 return a;
00212 }
00213
00214 a.n = 1;
00215 a.r = n;
00216 for (i=0; i<n; i++) a.mx[i] = 0.0;
00217 return a;
00218 }
00219
00220
00221
00231 imatrix make_imatrix1(int n)
00232 {
00233 int i;
00234 imatrix a;
00235
00236 a.n = a.r = 0;
00237 a.mx = NULL;
00238 a.sz = (int*)malloc(sizeof(int));
00239 if (a.sz==NULL) return a;
00240 a.sz[0] = n;
00241
00242 a.mx = (int*)malloc(n*sizeof(int));
00243 if (a.mx==NULL) {
00244 free(a.sz);
00245 return a;
00246 }
00247
00248 a.n = 1;
00249 a.r = n;
00250 for (i=0; i<n; i++) a.mx[i] = 0;
00251 return a;
00252 }
00253
00254
00255
00266 matrix make_matrix2(int n, int m)
00267 {
00268 int i, s;
00269 matrix a;
00270
00271 a.n = a.r = 0;
00272 a.mx = NULL;
00273 a.sz = (int*)malloc(2*sizeof(int));
00274 if (a.sz==NULL) return a;
00275 a.sz[0] = n;
00276 a.sz[1] = m;
00277 s = n*m;
00278
00279 a.mx = (double*)malloc(s*sizeof(double));
00280 if (a.mx==NULL) {
00281 free(a.sz);
00282 return a;
00283 }
00284
00285 a.n = 2;
00286 a.r = s;
00287 for (i=0; i<s; i++) a.mx[i] = 0.0;
00288 return a;
00289 }
00290
00291
00292
00303 imatrix make_imatrix2(int n, int m)
00304 {
00305 int i, s;
00306 imatrix a;
00307
00308 a.n = a.r = 0;
00309 a.mx = NULL;
00310 a.sz = (int*)malloc(2*sizeof(int));
00311 if (a.sz==NULL) return a;
00312 a.sz[0] = n;
00313 a.sz[1] = m;
00314 s = n*m;
00315
00316 a.mx = (int*)malloc(s*sizeof(int));
00317 if (a.mx==NULL) {
00318 free(a.sz);
00319 return a;
00320 }
00321
00322 a.n = 2;
00323 a.r = s;
00324 for (i=0; i<s; i++) a.mx[i] = 0;
00325 return a;
00326 }
00327
00328
00329
00340 matrix make_matrix(int n, int* sz)
00341 {
00342 int i, s;
00343 matrix a;
00344
00345 a.n = a.r = 0;
00346 a.sz = (int*)malloc(n*sizeof(int));
00347 if (a.sz==NULL) {
00348 a.mx = NULL;
00349 return a;
00350 }
00351 for (s=1,i=0; i<n; i++) {
00352 a.sz[i] = sz[i];
00353 s = s*sz[i];
00354 }
00355
00356 a.mx = (double*)malloc(s*sizeof(double));
00357 if (a.mx==NULL) {
00358 free(a.sz);
00359 return a;
00360 }
00361
00362 a.n = n;
00363 a.r = s;
00364 for (i=0; i<s; i++) a.mx[i] = 0.0;
00365 return a;
00366 }
00367
00368
00369
00380 imatrix make_imatrix(int n, int* sz)
00381 {
00382 int i, s;
00383 imatrix a;
00384
00385 a.n = a.r = 0;
00386 a.sz = (int*)malloc(n*sizeof(int));
00387 if (a.sz==NULL) {
00388 a.mx = NULL;
00389 return a;
00390 }
00391 for (s=1,i=0; i<n; i++) {
00392 a.sz[i] = sz[i];
00393 s = s*sz[i];
00394 }
00395
00396 a.mx = (int*)malloc(s*sizeof(int));
00397 if (a.mx==NULL) {
00398 free(a.sz);
00399 return a;
00400 }
00401
00402 a.n = n;
00403 a.r = s;
00404 for (i=0; i<s; i++) a.mx[i] = 0;
00405 return a;
00406 }
00407
00408
00409
00417 void free_matrix(matrix* a)
00418 {
00419 if(a->sz!=NULL) free(a->sz);
00420 if(a->mx!=NULL) free(a->mx);
00421 a->sz = NULL;
00422 a->mx = NULL;
00423 a->n = a->r = 0;
00424 }
00425
00426
00427
00435 void free_imatrix(imatrix* a)
00436 {
00437 if(a->sz!=NULL) free(a->sz);
00438 if(a->mx!=NULL) free(a->mx);
00439 a->sz = NULL;
00440 a->mx = NULL;
00441 a->n = a->r = 0;
00442 }
00443
00444
00445
00460 double* get_matrix(matrix mtx, ...)
00461 {
00462 int m, d;
00463 int* args;
00464 va_list argsptr;
00465
00466 if (mtx.n<1) return NULL;
00467 args = (int*)malloc(mtx.n*sizeof(int));
00468 if (args==NULL) return NULL;
00469
00470 va_start(argsptr, mtx);
00471 for (m=0; m<mtx.n; m++) {
00472 args[m] = (int)va_arg(argsptr, int);
00473 }
00474 va_end(argsptr);
00475
00476 int dx = args[0] - 1;
00477 for (d=1; d<mtx.n; d++) dx = dx*mtx.sz[d] + args[d] - 1;
00478 free(args);
00479
00480 if (dx>=mtx.r || dx<0) return NULL;
00481 return &mtx.mx[dx];
00482 }
00483
00484
00485
00500 int* get_imatrix(imatrix mtx, ...)
00501 {
00502 int m, d;
00503 int* args;
00504 va_list argsptr;
00505
00506 if (mtx.n<1) return NULL;
00507 args = (int*)malloc(mtx.n*sizeof(int));
00508 if (args==NULL) return NULL;
00509
00510 va_start(argsptr, mtx);
00511 for (m=0; m<mtx.n; m++) {
00512 args[m] = (int)va_arg(argsptr, int);
00513 }
00514 va_end(argsptr);
00515
00516 int dx = args[0] - 1;
00517 for (d=1; d<mtx.n; d++) dx = dx*mtx.sz[d] + args[d] - 1;
00518 free(args);
00519
00520 if (dx>=mtx.r || dx<0) return NULL;
00521 return &mtx.mx[dx];
00522 }
00523
00524
00525
00536 void copy_matrix(matrix src, matrix dst)
00537 {
00538 int i;
00539
00540 if (src.mx==NULL || dst.mx==NULL) return;
00541 if ((src.r!=dst.r)||(dst.n!=dst.n)) return;
00542
00543 for (i=0; i<src.n; i++) dst.sz[i] = src.sz[i];
00544 for (i=0; i<src.r; i++) dst.mx[i] = src.mx[i];
00545 }
00546
00547
00548
00559 void copy_imatrix(imatrix src, imatrix dst)
00560 {
00561 int i;
00562
00563 if (src.mx==NULL || dst.mx==NULL) return;
00564 if ((src.r!=dst.r)||(dst.n!=dst.n)) return;
00565
00566 for (i=0; i<src.n; i++) dst.sz[i] = src.sz[i];
00567 for (i=0; i<src.r; i++) dst.mx[i] = src.mx[i];
00568 }
00569
00570
00571
00582 matrix add_matrix(matrix a, matrix b)
00583 {
00584 int i;
00585 matrix c;
00586
00587 c.n = c.r = 0;
00588 c.sz = NULL;
00589 c.mx = NULL;
00590 if (a.mx==NULL || b.mx==NULL) return c;
00591 if (a.r != b.r) return c;
00592
00593 c = make_matrix(a.n, a.sz);
00594 for (i=0; i<c.r; i++) c.mx[i] = a.mx[i] + b.mx[i];
00595 return c;
00596 }
00597
00598
00599
00610 imatrix add_imatrix(imatrix a, imatrix b)
00611 {
00612 int i;
00613 imatrix c;
00614
00615 c.n = c.r = 0;
00616 c.sz = NULL;
00617 c.mx = NULL;
00618 if (a.mx==NULL || b.mx==NULL) return c;
00619 if (a.r != b.r) return c;
00620
00621 c = make_imatrix(a.n, a.sz);
00622 for (i=0; i<c.r; i++) c.mx[i] = a.mx[i] + b.mx[i];
00623 return c;
00624 }
00625
00626
00627
00638 matrix sub_matrix(matrix a, matrix b)
00639 {
00640 int i;
00641 matrix c;
00642
00643 c.n = c.r = 0;
00644 c.sz = NULL;
00645 c.mx = NULL;
00646 if (a.mx==NULL || b.mx==NULL) return c;
00647 if (a.r != b.r) return c;
00648
00649 c = make_matrix(a.n, a.sz);
00650 for (i=0; i<c.r; i++) c.mx[i] = a.mx[i] - b.mx[i];
00651 return c;
00652 }
00653
00654
00655
00666 imatrix sub_imatrix(imatrix a, imatrix b)
00667 {
00668 int i;
00669 imatrix c;
00670
00671 c.n = c.r = 0;
00672 c.sz = NULL;
00673 c.mx = NULL;
00674 if (a.mx==NULL || b.mx==NULL) return c;
00675 if (a.r != b.r) return c;
00676
00677 c = make_imatrix(a.n, a.sz);
00678 for (i=0; i<c.r; i++) c.mx[i] = a.mx[i] - b.mx[i];
00679 return c;
00680 }
00681
00682
00683
00694 matrix mlt_matrix(matrix a, matrix b)
00695 {
00696 int i, j, k, n, ii, aa, bb;
00697 int *sz, *sa, *sb, *sc, *cx;
00698 double st;
00699 matrix c;
00700
00701 c.n = c.r = 0;
00702 c.sz = NULL;
00703 c.mx = NULL;
00704 if (a.mx==NULL || b.mx==NULL) return c;
00705 if (a.sz[a.n-1]!=b.sz[0]) return c;
00706
00707 n = a.n + b.n - 2;
00708 sz = (int*)malloc(n*sizeof(int));
00709 if (sz==NULL) return c;
00710 sa = (int*)malloc(a.n*sizeof(int));
00711 if (sa==NULL) {free(sz); return c;}
00712 sb = (int*)malloc(b.n*sizeof(int));
00713 if (sb==NULL) {free(sz); free(sa); return c;}
00714 sc = (int*)malloc(n*sizeof(int));
00715 if (sc==NULL) {free(sz); free(sa); free(sb); return c;}
00716 cx = (int*)malloc(n*sizeof(int));
00717 if (cx==NULL) {free(sz); free(sa); free(sb); free(sc); return c;}
00718
00719 for (i=0; i<a.n-1; i++) sz[i] = a.sz[i];
00720 for (i=1; i<b.n; i++) sz[a.n-2+i] = b.sz[i];
00721
00722 sa[a.n-1] = sb[b.n-1] = sc[n-1] = 1;
00723 for (i=a.n-2; i>=0; i--) sa[i] = sa[i+1]*a.sz[i+1];
00724 for (i=b.n-2; i>=0; i--) sb[i] = sb[i+1]*b.sz[i+1];
00725 for (i=n-2; i>=0; i--) sc[i] = sc[i+1]*sz[i+1];
00726
00727 c = make_matrix(n, sz);
00728
00729 for (i=0; i<c.r; i++) {
00730 ii = i;
00731 for (j=0; j<c.n; j++) {
00732 cx[j] = ii / sc[j];
00733 ii = ii % sc[j];
00734 }
00735 aa = bb = 0;
00736 for (j=0; j<a.n-1; j++) aa = aa + sa[j]*cx[j];
00737 for (j=1; j<b.n; j++) bb = bb + sb[j]*cx[j+a.n-2];
00738
00739 st = 0.0;
00740 for (k=0; k<b.sz[0]; k++) st = st + a.mx[k+aa]*b.mx[bb+sb[0]*k];
00741 c.mx[i] = st;
00742 }
00743
00744 free(sz);
00745 free(sa);
00746 free(sb);
00747 free(sc);
00748 free(cx);
00749
00750 return c;
00751 }
00752
00753
00754
00765 imatrix mlt_imatrix(imatrix a, imatrix b)
00766 {
00767 int i, j, k, n, ii, aa, bb, st;
00768 int *sz, *sa, *sb, *sc, *cx;
00769 imatrix c;
00770
00771 c.n = c.r = 0;
00772 c.sz = NULL;
00773 c.mx = NULL;
00774 if (a.mx==NULL || b.mx==NULL) return c;
00775 if (a.sz[a.n-1]!=b.sz[0]) return c;
00776
00777 n = a.n + b.n - 2;
00778 sz = (int*)malloc(n*sizeof(int));
00779 if (sz==NULL) return c;
00780 sa = (int*)malloc(a.n*sizeof(int));
00781 if (sa==NULL) {free(sz); return c;}
00782 sb = (int*)malloc(b.n*sizeof(int));
00783 if (sb==NULL) {free(sz); free(sa); return c;}
00784 sc = (int*)malloc(n*sizeof(int));
00785 if (sc==NULL) {free(sz); free(sa); free(sb); return c;}
00786 cx = (int*)malloc(n*sizeof(int));
00787 if (cx==NULL) {free(sz); free(sa); free(sb); free(sc); return c;}
00788
00789
00790 for (i=0; i<a.n-1; i++) sz[i] = a.sz[i];
00791 for (i=1; i<b.n; i++) sz[a.n-2+i] = b.sz[i];
00792
00793 sa[a.n-1] = sb[b.n-1] = sc[n-1] = 1;
00794 for (i=a.n-2; i>=0; i--) sa[i] = sa[i+1]*a.sz[i+1];
00795 for (i=b.n-2; i>=0; i--) sb[i] = sb[i+1]*b.sz[i+1];
00796 for (i=n-2; i>=0; i--) sc[i] = sc[i+1]*sz[i+1];
00797
00798 c = make_imatrix(n, sz);
00799
00800 for (i=0; i<c.r; i++) {
00801 ii = i;
00802 for (j=0; j<c.n; j++) {
00803 cx[j] = ii / sc[j];
00804 ii = ii % sc[j];
00805 }
00806 aa = bb = 0;
00807 for (j=0; j<a.n-1; j++) aa = aa + sa[j]*cx[j];
00808 for (j=1; j<b.n; j++) bb = bb + sb[j]*cx[j+a.n-2];
00809
00810 st = 0;
00811 for (k=0; k<b.sz[0]; k++) st = st + a.mx[k+aa]*b.mx[bb+sb[0]*k];
00812 c.mx[i] = st;
00813 }
00814
00815 free(sz);
00816 free(sa);
00817 free(sb);
00818 free(sc);
00819 free(cx);
00820
00821 return c;
00822 }
00823
00824
00825
00834 void print_matrix(FILE* fp, matrix a)
00835 {
00836 int i;
00837 for (i=0; i<a.r; i++) {
00838 fprintf(fp, " %15lf", a.mx[i]);
00839 if ((i+1)%a.sz[a.n-1]==0) fprintf(fp, "\n");
00840 }
00841
00842 }
00843
00844
00845
00854 void print_imatrix(FILE* fp, imatrix a)
00855 {
00856 int i;
00857 for (i=0; i<a.r; i++) {
00858 fprintf(fp, " %6d", a.mx[i]);
00859 if ((i+1)%a.sz[a.n-1]==0) fprintf(fp, "\n");
00860 }
00861 fprintf(stdout,"\n");
00862 }
00863
00864
00865
00881 matrix decompQR(matrix xx, imatrix col)
00882 {
00883 int i, j, r, n, m, sz[2];
00884 double s, t, u;
00885 matrix nsq, isq, x, a;
00886
00887 a.n = a.r = 0;
00888 a.sz = NULL;
00889 a.mx = NULL;
00890 if (xx.mx==NULL || col.mx==NULL) return a;
00891 if (xx.n!=2 || (xx.sz[1]!=col.sz[0])) return a;
00892
00893 nsq.mx = isq.mx = x.mx = NULL;
00894 n = xx.sz[0];
00895 m = xx.sz[1];
00896 sz[0] = sz[1] = m;
00897 nsq = make_matrix(1, &m);
00898 isq = make_matrix(1, &m);
00899 x = make_matrix(xx.n, xx.sz);
00900 a = make_matrix(xx.n, sz);
00901 if (nsq.mx==NULL || isq.mx==NULL || x.mx==NULL || a.mx==NULL) {
00902 free_matrix(&nsq);
00903 free_matrix(&isq);
00904 free_matrix(&x);
00905 return a;
00906 }
00907
00908 for (i=0; i<xx.r; i++) x.mx[i] = xx.mx[i];
00909
00910 for (i=1; i<=m; i++) {
00911 for (s=0.0,j=1; j<=n; j++) s = s + Mx(x, j, i)*Mx(x, j, i);
00912 Vt(nsq, i) = s;
00913 Vt(isq, i) = ((Vt(nsq, i)!=0) ? Vt(nsq,i) : -1.0);
00914 Vt(col, i) = i;
00915 }
00916
00917 for (r=1; r<=m; r++) {
00918 if (r!=1) {
00919 j = r; u = 0.0;
00920 for (i=r; i<=m; i++) {
00921 t = Vt(nsq,i)/Vt(isq,i);
00922 if (t>u) { u = t; j = i; }
00923 }
00924 i = Vt(col,j); Vt(col,j) = Vt(col,r); Vt(col,r) = i;
00925 t = Vt(nsq,j); Vt(nsq,j) = Vt(nsq,r); Vt(nsq,r) = t;
00926 t = Vt(isq,j); Vt(isq,j) = Vt(isq,r); Vt(isq,r) = t;
00927 for (i=1; i<=n; i++) {
00928 t = Mx(x,i,j);
00929 Mx(x,i,j) = Mx(x,i,r);
00930 Mx(x,i,r) = t;
00931 }
00932 }
00933
00934 for (u=0.0,i=r; i<=n; i++) u = u + Mx(x,i,r)*Mx(x,i,r);
00935 u = sqrt(u);
00936 if (Mx(x,r,r)<0.0) u = -u;
00937 Mx(x, r, r) = Mx(x,r,r) + u;
00938 t = 1.0/(Mx(x,r,r)*u);
00939
00940 for (j=1; j<=r-1; j++) Mx(x,r,j)=0.0;
00941 for (j=r+1; j<=m; j++) {
00942 for (s=0.0,i=r; i<=n; i++) s = s + Mx(x,i,r)*Mx(x,i,j);
00943 for (i=r; i<=n; i++) Mx(x,i,j) = Mx(x,i,j) - s*t*Mx(x,i,r);
00944 }
00945 Mx(x,r,r) = -u;
00946 }
00947
00948 for (i=1; i<=m; i++) {
00949 for (j=1; j<=m; j++) Mx(a,i,j) = Mx(x,i,j);
00950 }
00951
00952 free_matrix(&nsq);
00953 free_matrix(&isq);
00954 free_matrix(&x);
00955
00956 return a;
00957 }
00958
00959
00960
00980 matrix minimum2(matrix y, matrix x)
00981 {
00982 int i, m;
00983 imatrix cl;
00984 matrix rx, rt, rr, aa, bb, cc;
00985
00986 cc.n = cc.r = 0;
00987 cc.sz = NULL;
00988 cc.mx = NULL;
00989 if (y.mx==NULL || x.mx==NULL) return cc;
00990 if (y.sz[0]!=x.sz[0]) return cc;
00991
00992
00993 m = x.sz[1];
00994 cl = make_imatrix(1, &m);
00995 if (cl.mx==NULL) return cc;
00996 cc = make_matrix (1, &m);
00997 if (cc.mx==NULL) {free_imatrix(&cl); return cc;}
00998
00999 rx = decompQR(x, cl);
01000 rr = invrU_matrix(rx),
01001 rt = mlt_matrix(rr, trans_matrix(rr));
01002 aa = mlt_matrix(trans_matrix(x), y),
01003 bb = mlt_matrix(rt, aa);
01004
01005 if (bb.mx!=NULL) for (i=0; i<m; i++) cc.mx[cl.mx[i]-1] = bb.mx[i];
01006
01007 free_imatrix(&cl);
01008 free_matrix(&rx);
01009 free_matrix(&rr);
01010 free_matrix(&rt);
01011 free_matrix(&aa);
01012 free_matrix(&bb);
01013
01014 return cc;
01015 }
01016
01017
01018
01027 matrix trans_matrix(matrix a)
01028 {
01029 int i, j, k;
01030 matrix c;
01031 int sz[2];
01032
01033 c.n = c.r = 0;
01034 c.sz = NULL;
01035 c.mx = NULL;
01036 if (a.mx==NULL || a.n!=2) return c;
01037
01038 sz[0] = a.sz[1];
01039 sz[1] = a.sz[0];
01040 c = make_matrix(a.n, sz);
01041
01042 for (k=0; k<c.r; k++) {
01043 i = k/a.sz[1];
01044 j = k%a.sz[1];
01045 c.mx[j*sz[1]+i] = a.mx[k];
01046 }
01047 return c;
01048 }
01049
01050
01051
01060 matrix invrU_matrix(matrix x)
01061 {
01062 int i, j, k, n;
01063 double t, u, det;
01064 matrix a;
01065
01066 a.n = a.r = 0;
01067 a.sz = NULL;
01068 a.mx = NULL;
01069 if (x.mx==NULL) return a;
01070 if (x.sz[0]!=x.sz[1]) return a;
01071
01072 n = x.sz[0];
01073 for (j=1; j<n; j++) {
01074 for (i=j+1; i<=n; i++) {
01075 if (Mx(x,i,j) != 0.0) return a;
01076 }
01077 }
01078
01079 det = 1.0;
01080 a = make_matrix(x.n, x.sz);
01081 if (a.mx==NULL) return a;
01082
01083 for (i=0; i<a.r; i++) a.mx[i] = x.mx[i];
01084
01085 for (k=1; k<=n; k++) {
01086 t = Mx(a,k,k);
01087 det = det*t;
01088 for (i=1; i<=n; i++) Mx(a,i,k) = Mx(a,i,k)/t;
01089 Mx(a,k,k) = 1.0/t;
01090
01091 for (j=1; j<=n; j++) {
01092 if (j!=k) {
01093 u = Mx(a,k,j);
01094 for (i=1; i<=n; i++) {
01095 if (i!=k) Mx(a,i,j) = Mx(a,i,j) - Mx(a,i,k)*u;
01096 else Mx(a,i,j) = - u/t;
01097 }
01098 }
01099 }
01100 }
01101 return a;
01102 }
01103
