Embedded Files
> Back
/* Library for Matrix Operations --- Extended */
#include "matrix0.c"
mat_t *mat_copy(mat_t *org, mat_t *mat) { int r = org->rows, c = org->cols, ctr = r * c; double *o = org->elem, *m = mat_size(mat, r, c)->elem; for(; ctr > 0; ctr--) *m++ = *o++; return mat;}
mat_t *mat_trans(mat_t *org, mat_t *mat) { int r = org->rows, c = org->cols, i = 1, j; double *o = org->elem, *t_base = mat_size(mat, c, r)->elem, *t; for(; i <= r; i++, t_base++) for(t = t_base, j = 1; j <= c; j++, t += r) *t = *o++; return mat;}
mat_t *mat_sprod(mat_t *org, double coef, mat_t *mat) { int r = org->rows, c = org->cols, ctr = r * c; double *o = org->elem, *m = mat_size(mat, r, c)->elem; for(; ctr > 0; ctr--) *m++ = *o++ * coef; return mat;}
mat_t *mat_add(mat_t *frst, mat_t *scnd, mat_t *mat) { int r = frst->rows, c = frst->cols, scnd_r = scnd->rows, scnd_c = scnd->cols, ctr; double *f, *s, *m; if((r != scnd_r) || (c != scnd_c)) { printf("*** DIM MISMATCH ***"); exit(1); } for(f = frst->elem, s = scnd->elem, m = mat_size(mat, r, c)->elem, ctr = r * c; ctr > 0; ctr--) *m++ = *f++ + *s++; return mat;}
mat_t *mat_sub(mat_t *frst, mat_t *scnd, mat_t *mat) { int r = frst->rows, c = frst->cols, scnd_r = scnd->rows, scnd_c = scnd->cols, ctr; double *f, *s, *m; if((r != scnd_r) || (c != scnd_c)) { printf("*** DIM MISMATCH ***"); exit(1); } for(f = frst->elem, s = scnd->elem, m = mat_size(mat, r, c)->elem, ctr = r * c; ctr > 0; ctr--) *m++ = *f++ - *s++; return mat;}
mat_t *mat_mul(mat_t *frst, mat_t *scnd, mat_t *mat) { int r, c, frst_c = frst->cols, scnd_r = scnd->rows, i = 1, j, ctr; double sum, *f_base, *s_base, *f, *s, *m; if(frst_c != scnd_r) { printf("*** DIM MISMATCH ***"); exit(1); } m = mat_size(mat, r = frst->rows, c = scnd->cols)->elem; for(f_base = frst->elem; i <= r; i++, f_base += frst_c) { for(s_base = scnd->elem, j = 1; j <= c; j++, s_base++) { sum = 0; for(f = f_base, s = s_base, ctr = frst_c; ctr > 0; ctr--, s += c) sum += *f++ * *s; *m++ = sum; } } return mat;}
mat_t *mat_minor(mat_t *org, int row, int col, mat_t *mat) { int r = org->rows, c = org->cols, i = 1, j; double *o = org->elem, *m = mat_size(mat, r - 1, c - 1)->elem; for(; i <= r; i++) { if(i != row) { for(j = 1; j <= c; j++, o++) if(j != col) *m++ = *o; } else o += c; } return mat;}
double mat_det(mat_t *mat) { mat_t tmp; int r = mat->rows, c = mat->cols, j = 1, sgn = 1; double sum = 0, *m = mat->elem; if(r != c) { printf("*** NOT SQUARE ***"); exit(1); } if(r == 1) return *(mat->elem); mat_init(&tmp); for(; j <= c; j++, m++) { if(sgn) { if(*m != 0) sum += *m * mat_det(mat_minor(mat, 1, j, &tmp)); sgn = 0; } else { if(*m != 0) sum -= *m * mat_det(mat_minor(mat, 1, j, &tmp)); sgn = 1; } /* 要素が 0 のときのみ縮小行列の行列式を再帰的に算出 */ } mat_free(&tmp); return sum;}
double mat_adjelem(mat_t *mat, int row, int col) { mat_t tmp; double det = mat_det(mat_minor(mat, row, col, mat_init(&tmp))); mat_free(&tmp); return ((row + col) % 2) ? -det : det;}
mat_t *mat_adj(mat_t *org, mat_t *mat) { mat_t tmp; int r = org->rows, c = org->cols, i = 1, j; double *t = mat_size(mat_init(&tmp), r, c)->elem; for(; i <= r; i++) for(j = 1; j <= c; j++) *t++ = mat_adjelem(org, i, j); mat_trans(&tmp, mat); mat_free(&tmp); return mat;}
double mat_inv(mat_t *org, mat_t *mat) { mat_t tmp; double det = mat_det(org); if(nearint(det) == 0) return 0; if(org->rows == 1) { mat_write(mat_size(mat, 1, 1), 1, 1, (double) 1 / det); return det; } mat_sprod(mat_adj(org, mat_init(&tmp)), (double) 1 / det, mat_size(mat, org->rows, org->cols)); mat_free(&tmp); return det;}
double mat_eq(mat_t *coef, mat_t *cons, mat_t *mat) { mat_t tmp; double det = mat_inv(coef, mat_init(&tmp)); if(det == 0) return det; mat_mul(&tmp, cons, mat_size(mat, coef->rows, coef->cols + 1)); mat_free(&tmp); return det;}
int mat_govrow(mat_t *mat, int col, double *govcoef) { int r = mat->rows, c, govrow = col, i = col; double abscoef, tmp, abstmp, *m; if((col < 1) || (r < col)) { printf("*** OUT OF DIM ***"); exit(1); } m = mat->elem + ((c = mat->cols) + 1) * (i - 1); abscoef = fabs(*govcoef = *m); for(i++, m += c; i <= r; i++, m += c) { if((abstmp = fabs(tmp = *m)) > abscoef) { govrow = i; *govcoef = tmp; abscoef = abstmp; } } return govrow;}
#include "matrix0.c"
mat_t *mat_copy(mat_t *org, mat_t *mat) { int r = org->rows, c = org->cols, ctr = r * c; double *o = org->elem, *m = mat_size(mat, r, c)->elem; for(; ctr > 0; ctr--) *m++ = *o++; return mat;}
mat_t *mat_trans(mat_t *org, mat_t *mat) { int r = org->rows, c = org->cols, i = 1, j; double *o = org->elem, *t_base = mat_size(mat, c, r)->elem, *t; for(; i <= r; i++, t_base++) for(t = t_base, j = 1; j <= c; j++, t += r) *t = *o++; return mat;}
mat_t *mat_sprod(mat_t *org, double coef, mat_t *mat) { int r = org->rows, c = org->cols, ctr = r * c; double *o = org->elem, *m = mat_size(mat, r, c)->elem; for(; ctr > 0; ctr--) *m++ = *o++ * coef; return mat;}
mat_t *mat_add(mat_t *frst, mat_t *scnd, mat_t *mat) { int r = frst->rows, c = frst->cols, scnd_r = scnd->rows, scnd_c = scnd->cols, ctr; double *f, *s, *m; if((r != scnd_r) || (c != scnd_c)) { printf("*** DIM MISMATCH ***"); exit(1); } for(f = frst->elem, s = scnd->elem, m = mat_size(mat, r, c)->elem, ctr = r * c; ctr > 0; ctr--) *m++ = *f++ + *s++; return mat;}
mat_t *mat_sub(mat_t *frst, mat_t *scnd, mat_t *mat) { int r = frst->rows, c = frst->cols, scnd_r = scnd->rows, scnd_c = scnd->cols, ctr; double *f, *s, *m; if((r != scnd_r) || (c != scnd_c)) { printf("*** DIM MISMATCH ***"); exit(1); } for(f = frst->elem, s = scnd->elem, m = mat_size(mat, r, c)->elem, ctr = r * c; ctr > 0; ctr--) *m++ = *f++ - *s++; return mat;}
mat_t *mat_mul(mat_t *frst, mat_t *scnd, mat_t *mat) { int r, c, frst_c = frst->cols, scnd_r = scnd->rows, i = 1, j, ctr; double sum, *f_base, *s_base, *f, *s, *m; if(frst_c != scnd_r) { printf("*** DIM MISMATCH ***"); exit(1); } m = mat_size(mat, r = frst->rows, c = scnd->cols)->elem; for(f_base = frst->elem; i <= r; i++, f_base += frst_c) { for(s_base = scnd->elem, j = 1; j <= c; j++, s_base++) { sum = 0; for(f = f_base, s = s_base, ctr = frst_c; ctr > 0; ctr--, s += c) sum += *f++ * *s; *m++ = sum; } } return mat;}
mat_t *mat_minor(mat_t *org, int row, int col, mat_t *mat) { int r = org->rows, c = org->cols, i = 1, j; double *o = org->elem, *m = mat_size(mat, r - 1, c - 1)->elem; for(; i <= r; i++) { if(i != row) { for(j = 1; j <= c; j++, o++) if(j != col) *m++ = *o; } else o += c; } return mat;}
double mat_det(mat_t *mat) { mat_t tmp; int r = mat->rows, c = mat->cols, j = 1, sgn = 1; double sum = 0, *m = mat->elem; if(r != c) { printf("*** NOT SQUARE ***"); exit(1); } if(r == 1) return *(mat->elem); mat_init(&tmp); for(; j <= c; j++, m++) { if(sgn) { if(*m != 0) sum += *m * mat_det(mat_minor(mat, 1, j, &tmp)); sgn = 0; } else { if(*m != 0) sum -= *m * mat_det(mat_minor(mat, 1, j, &tmp)); sgn = 1; } /* 要素が 0 のときのみ縮小行列の行列式を再帰的に算出 */ } mat_free(&tmp); return sum;}
double mat_adjelem(mat_t *mat, int row, int col) { mat_t tmp; double det = mat_det(mat_minor(mat, row, col, mat_init(&tmp))); mat_free(&tmp); return ((row + col) % 2) ? -det : det;}
mat_t *mat_adj(mat_t *org, mat_t *mat) { mat_t tmp; int r = org->rows, c = org->cols, i = 1, j; double *t = mat_size(mat_init(&tmp), r, c)->elem; for(; i <= r; i++) for(j = 1; j <= c; j++) *t++ = mat_adjelem(org, i, j); mat_trans(&tmp, mat); mat_free(&tmp); return mat;}
double mat_inv(mat_t *org, mat_t *mat) { mat_t tmp; double det = mat_det(org); if(nearint(det) == 0) return 0; if(org->rows == 1) { mat_write(mat_size(mat, 1, 1), 1, 1, (double) 1 / det); return det; } mat_sprod(mat_adj(org, mat_init(&tmp)), (double) 1 / det, mat_size(mat, org->rows, org->cols)); mat_free(&tmp); return det;}
double mat_eq(mat_t *coef, mat_t *cons, mat_t *mat) { mat_t tmp; double det = mat_inv(coef, mat_init(&tmp)); if(det == 0) return det; mat_mul(&tmp, cons, mat_size(mat, coef->rows, coef->cols + 1)); mat_free(&tmp); return det;}
int mat_govrow(mat_t *mat, int col, double *govcoef) { int r = mat->rows, c, govrow = col, i = col; double abscoef, tmp, abstmp, *m; if((col < 1) || (r < col)) { printf("*** OUT OF DIM ***"); exit(1); } m = mat->elem + ((c = mat->cols) + 1) * (i - 1); abscoef = fabs(*govcoef = *m); for(i++, m += c; i <= r; i++, m += c) { if((abstmp = fabs(tmp = *m)) > abscoef) { govrow = i; *govcoef = tmp; abscoef = abstmp; } } return govrow;}
Google Sites
Report abuse