![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mat.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to mat.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/mat.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3:
! 4: void addmat(vl,a,b,c)
! 5: VL vl;
! 6: MAT a,b,*c;
! 7: {
! 8: int row,col,i,j;
! 9: MAT t;
! 10: pointer *ab,*bb,*tb;
! 11:
! 12: if ( !a )
! 13: *c = b;
! 14: else if ( !b )
! 15: *c = a;
! 16: else if ( (a->row != b->row) || (a->col != b->col) ) {
! 17: *c = 0; error("addmat : size mismatch");
! 18: } else {
! 19: row = a->row; col = a->col;
! 20: MKMAT(t,row,col);
! 21: for ( i = 0; i < row; i++ )
! 22: for ( j = 0, ab = BDY(a)[i], bb = BDY(b)[i], tb = BDY(t)[i];
! 23: j < col; j++ )
! 24: addr(vl,(Obj)ab[j],(Obj)bb[j],(Obj *)&tb[j]);
! 25: *c = t;
! 26: }
! 27: }
! 28:
! 29: void submat(vl,a,b,c)
! 30: VL vl;
! 31: MAT a,b,*c;
! 32: {
! 33: int row,col,i,j;
! 34: MAT t;
! 35: pointer *ab,*bb,*tb;
! 36:
! 37: if ( !a )
! 38: chsgnmat(b,c);
! 39: else if ( !b )
! 40: *c = a;
! 41: else if ( (a->row != b->row) || (a->col != b->col) ) {
! 42: *c = 0; error("submat : size mismatch");
! 43: } else {
! 44: row = a->row; col = a->col;
! 45: MKMAT(t,row,col);
! 46: for ( i = 0; i < row; i++ )
! 47: for ( j = 0, ab = BDY(a)[i], bb = BDY(b)[i], tb = BDY(t)[i];
! 48: j < col; j++ )
! 49: subr(vl,(Obj)ab[j],(Obj)bb[j],(Obj *)&tb[j]);
! 50: *c = t;
! 51: }
! 52: }
! 53:
! 54: void mulmat(vl,a,b,c)
! 55: VL vl;
! 56: Obj a,b,*c;
! 57: {
! 58: if ( !a || !b )
! 59: *c = 0;
! 60: else if ( OID(a) <= O_R )
! 61: mulrmat(vl,(Obj)a,(MAT)b,(MAT *)c);
! 62: else if ( OID(b) <= O_R )
! 63: mulrmat(vl,(Obj)b,(MAT)a,(MAT *)c);
! 64: else
! 65: switch ( OID(a) ) {
! 66: case O_VECT:
! 67: switch ( OID(b) ) {
! 68: case O_MAT:
! 69: mulvectmat(vl,(VECT)a,(MAT)b,(VECT *)c); break;
! 70: case O_VECT: default:
! 71: notdef(vl,a,b,c); break;
! 72: }
! 73: break;
! 74: case O_MAT:
! 75: switch ( OID(b) ) {
! 76: case O_VECT:
! 77: mulmatvect(vl,(MAT)a,(VECT)b,(VECT *)c); break;
! 78: case O_MAT:
! 79: mulmatmat(vl,(MAT)a,(MAT)b,(MAT *)c); break;
! 80: default:
! 81: notdef(vl,a,b,c); break;
! 82: }
! 83: break;
! 84: default:
! 85: notdef(vl,a,b,c); break;
! 86: }
! 87: }
! 88:
! 89: void divmat(vl,a,b,c)
! 90: VL vl;
! 91: Obj a,b,*c;
! 92: {
! 93: Obj t;
! 94:
! 95: if ( !b )
! 96: error("divmat : division by 0");
! 97: else if ( !a )
! 98: *c = 0;
! 99: else if ( OID(b) > O_R )
! 100: notdef(vl,a,b,c);
! 101: else {
! 102: divr(vl,(Obj)ONE,b,&t); mulrmat(vl,t,(MAT)a,(MAT *)c);
! 103: }
! 104: }
! 105:
! 106: void chsgnmat(a,b)
! 107: MAT a,*b;
! 108: {
! 109: MAT t;
! 110: int row,col,i,j;
! 111: pointer *ab,*tb;
! 112:
! 113: if ( !a )
! 114: *b = 0;
! 115: else {
! 116: row = a->row; col = a->col;
! 117: MKMAT(t,row,col);
! 118: for ( i = 0; i < row; i++ )
! 119: for ( j = 0, ab = BDY(a)[i], tb = BDY(t)[i];
! 120: j < col; j++ )
! 121: chsgnr((Obj)ab[j],(Obj *)&tb[j]);
! 122: *b = t;
! 123: }
! 124: }
! 125:
! 126: void pwrmat(vl,a,r,c)
! 127: VL vl;
! 128: MAT a;
! 129: Obj r;
! 130: MAT *c;
! 131: {
! 132: if ( !a )
! 133: *c = 0;
! 134: else if ( !r || !NUM(r) || !RATN(r) ||
! 135: !INT(r) || (SGN((Q)r)<0) || (PL(NM((Q)r))>1) ) {
! 136: *c = 0; error("pwrmat : invalid exponent");
! 137: } else if ( a->row != a->col ) {
! 138: *c = 0; error("pwrmat : non square matrix");
! 139: } else
! 140: pwrmatmain(vl,a,QTOS((Q)r),c);
! 141: }
! 142:
! 143: void pwrmatmain(vl,a,e,c)
! 144: VL vl;
! 145: MAT a;
! 146: int e;
! 147: MAT *c;
! 148: {
! 149: MAT t,s;
! 150:
! 151: if ( e == 1 ) {
! 152: *c = a;
! 153: return;
! 154: }
! 155:
! 156: pwrmatmain(vl,a,e/2,&t);
! 157: mulmat(vl,(Obj)t,(Obj)t,(Obj *)&s);
! 158: if ( e % 2 )
! 159: mulmat(vl,(Obj)s,(Obj)a,(Obj *)c);
! 160: else
! 161: *c = s;
! 162: }
! 163:
! 164: void mulrmat(vl,a,b,c)
! 165: VL vl;
! 166: Obj a;
! 167: MAT b,*c;
! 168: {
! 169: int row,col,i,j;
! 170: MAT t;
! 171: pointer *bb,*tb;
! 172:
! 173: if ( !a || !b )
! 174: *c = 0;
! 175: else {
! 176: row = b->row; col = b->col;
! 177: MKMAT(t,row,col);
! 178: for ( i = 0; i < row; i++ )
! 179: for ( j = 0, bb = BDY(b)[i], tb = BDY(t)[i];
! 180: j < col; j++ )
! 181: mulr(vl,(Obj)a,(Obj)bb[j],(Obj *)&tb[j]);
! 182: *c = t;
! 183: }
! 184: }
! 185:
! 186: void mulmatmat(vl,a,b,c)
! 187: VL vl;
! 188: MAT a,b,*c;
! 189: {
! 190: int arow,bcol,i,j,k,m;
! 191: MAT t;
! 192: pointer s,u,v;
! 193: pointer *ab,*tb;
! 194:
! 195: if ( a->col != b->row ) {
! 196: *c = 0; error("mulmat : size mismatch");
! 197: } else {
! 198: arow = a->row; m = a->col; bcol = b->col;
! 199: MKMAT(t,arow,bcol);
! 200: for ( i = 0; i < arow; i++ )
! 201: for ( j = 0, ab = BDY(a)[i], tb = BDY(t)[i]; j < bcol; j++ ) {
! 202: for ( k = 0, s = 0; k < m; k++ ) {
! 203: mulr(vl,(Obj)ab[k],(Obj)BDY(b)[k][j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
! 204: }
! 205: tb[j] = s;
! 206: }
! 207: *c = t;
! 208: }
! 209: }
! 210:
! 211: void mulmatvect(vl,a,b,c)
! 212: VL vl;
! 213: MAT a;
! 214: VECT b;
! 215: VECT *c;
! 216: {
! 217: int arow,i,j,m;
! 218: VECT t;
! 219: pointer s,u,v;
! 220: pointer *ab;
! 221:
! 222: if ( !a || !b )
! 223: *c = 0;
! 224: else if ( a->col != b->len ) {
! 225: *c = 0; error("mulmatvect : size mismatch");
! 226: } else {
! 227: arow = a->row; m = a->col;
! 228: MKVECT(t,arow);
! 229: for ( i = 0; i < arow; i++ ) {
! 230: for ( j = 0, s = 0, ab = BDY(a)[i]; j < m; j++ ) {
! 231: mulr(vl,(Obj)ab[j],(Obj)BDY(b)[j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
! 232: }
! 233: BDY(t)[i] = s;
! 234: }
! 235: *c = t;
! 236: }
! 237: }
! 238:
! 239: void mulvectmat(vl,a,b,c)
! 240: VL vl;
! 241: VECT a;
! 242: MAT b;
! 243: VECT *c;
! 244: {
! 245: int bcol,i,j,m;
! 246: VECT t;
! 247: pointer s,u,v;
! 248:
! 249: if ( !a || !b )
! 250: *c = 0;
! 251: else if ( a->len != b->row ) {
! 252: *c = 0; error("mulvectmat : size mismatch");
! 253: } else {
! 254: bcol = b->col; m = a->len;
! 255: MKVECT(t,bcol);
! 256: for ( j = 0; j < bcol; j++ ) {
! 257: for ( i = 0, s = 0; i < m; i++ ) {
! 258: mulr(vl,(Obj)BDY(a)[i],(Obj)BDY(b)[i][j],(Obj *)&u); addr(vl,(Obj)s,(Obj)u,(Obj *)&v); s = v;
! 259: }
! 260: BDY(t)[j] = s;
! 261: }
! 262: *c = t;
! 263: }
! 264: }
! 265:
! 266: int compmat(vl,a,b)
! 267: VL vl;
! 268: MAT a,b;
! 269: {
! 270: int i,j,t,row,col;
! 271:
! 272: if ( !a )
! 273: return b?-1:0;
! 274: else if ( !b )
! 275: return 1;
! 276: else if ( a->row != b->row )
! 277: return a->row>b->row ? 1 : -1;
! 278: else if (a->col != b->col )
! 279: return a->col > b->col ? 1 : -1;
! 280: else {
! 281: row = a->row; col = a->col;
! 282: for ( i = 0; i < row; i++ )
! 283: for ( j = 0; j < col; j++ )
! 284: if ( t = compr(vl,(Obj)BDY(a)[i][j],(Obj)BDY(b)[i][j]) )
! 285: return t;
! 286: return 0;
! 287: }
! 288: }
! 289:
! 290: pointer **almat_pointer(n,m)
! 291: int n,m;
! 292: {
! 293: pointer **mat;
! 294: int i;
! 295:
! 296: mat = (pointer **)MALLOC(n*sizeof(pointer *));
! 297: for ( i = 0; i < n; i++ )
! 298: mat[i] = (pointer *)CALLOC(m,sizeof(pointer));
! 299: return mat;
! 300: }