Mr Nobody
BANNED | Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору #include <stdio.h> #include <alloc.h> #include <math.h> int cholesky_decomp(float *mat_a[], float *mat_b[], int sz) { float sum; int ii, jj, kk; for(ii = 0;ii < sz;ii++) { for(jj = ii;jj < sz;jj++) { sum = mat_a[ii][jj]; for(kk = ii - 1;kk >= 0;kk--) sum = sum - mat_a[ii][kk]*mat_a[jj][kk]; if(ii == jj) { if(sum <= .0F) { printf("Cholesky's Method failed!\n"); return 0; } mat_b[ii][ii] = sqrt(sum); } else mat_b[jj][ii] = mat_a[jj][ii] = sum/mat_b[ii][ii]; } } return 1; } void cholesky_solver(float *mat_b[], float right[], float xx[], int sz) { float *p_temp, sum; int ii, jj; p_temp = (float*)malloc(sz*sizeof(right[0])); for(ii = 0; ii < 4; ii++) { sum = 0.0f; for(jj = 0; jj < ii; jj++) sum += mat_b[ii][jj]*p_temp[jj]; p_temp[ii] = (right[ii] - sum)/mat_b[ii][ii]; } for(ii = 3; ii >= 0; ii--) { sum = 0.0f; for(jj = 3; jj > ii; jj--) sum += mat_b[jj][ii]*xx[jj]; xx[ii] = (p_temp[ii] - sum)/mat_b[ii][ii]; } free(p_temp); } float left[4][4] = {{1.0f, 1.0f, 1.0f, 1.0f}, // This matrix generate by pascal(4) from Matlab {1.0f, 2.0f, 3.0f, 4.0f}, {1.0f, 3.0f, 6.0f, 10.0f}, {1.0f, 4.0f, 10.0f, 20.0f}}; float bb[4][4]; float *p_left[4] = {(float*)&left[0], (float*)&left[1], (float*)&left[2], (float*)&left[3]}; float *p_bb[4] = {(float*)&bb[0], (float*)&bb[1], (float*)&bb[2], (float*)&bb[3]}; float right[4] = {6.0f, 16.0f, 33.0f, 59.0f}; float xx[4]; int main() { int ii; if(!cholesky_decomp(p_left, p_bb, 4)) return 0; cholesky_solver(p_bb, right, xx, 4); printf("Cholesky's Method\n\n"); for(ii = 0; ii < 4; ii++) printf(" X%d = %f\n", ii, xx[ii]); return 0; } | Всего записей: 350 | Зарегистр. 19-09-2007 | Отправлено: 12:19 31-10-2007 | Исправлено: Mr Nobody, 12:25 31-10-2007 |
|