/* Lab. Calcolo II - Esempio di codice Example code: Det calculation with pivoting. (20000814 prelz@mi.infn.it) */ #include <stdio.h> #include <stdlib.h> #include <math.h> #define MATRIX_SIZE 6 #define RANDOM_SEED 1515 #define N_REPEATS 50 #ifndef TRUE #define TRUE 1 #endif #ifndef FALSE #define FALSE 0 #endif double determinant(double **matrix, int n, int do_pivot); int main(int argc, char *argv[]) { double matrix[MATRIX_SIZE][MATRIX_SIZE]; double *m_pass[MATRIX_SIZE]; double result_piv, result_nopiv; int i,j,k; /* Initialize m_pass: it is used to pass matrix as a double** */ for (i=0;i<MATRIX_SIZE;i++) m_pass[i]=matrix[i]; srand(RANDOM_SEED); for (k=0; k<N_REPEATS; k++) { for (i=0;i<MATRIX_SIZE;i++) for (j=0;j<MATRIX_SIZE;j++) { matrix[i][j] = (double)(rand()%10000000); if (matrix[i][j] < 5000000) matrix[i][j] = 1/matrix[i][j]; } result_nopiv = determinant(m_pass, MATRIX_SIZE, FALSE); result_piv = determinant(m_pass, MATRIX_SIZE, TRUE); if (fabs(result_nopiv - result_piv) > 1E-3) { print_matrix(m_pass, MATRIX_SIZE); printf ("Determinant without pivoting is %30.25g\n",result_nopiv); printf ("Determinant with pivoting is %30.25g\n\n",result_piv); } } exit(0); } double determinant(double **matrix, int n, int do_pivot) { int i,j,k; int max_row, max_column; double **m; /* Temporary copy that we can manipulate */ double *swap; double result; double sign = 1; double largest; double mul_factor; swap = (double *)malloc(n * sizeof(double)); if (swap == NULL) return (-1); m = (double **)malloc(n * sizeof(double *)); if (m == NULL) { free(swap); return(-1); } for (i=0;i<=n;i++) { m[i] = (double *)malloc(n * sizeof(double)); if (m[i] == NULL) { free(swap); for(j=0;j<i;j++) free(m[j]); free(m); return(-1); } memcpy(m[i],matrix[i],n * sizeof(double)); } for (i=0;i<n;i++) { if (do_pivot) { /* Find a good pivot element in the lower right (n-i)x(n-i) */ /* minor to put in the upper left corner. */ /* We'll choose the largest element in the minor. */ largest = 0; for (j=i; j<n; j++) { for (k=i; k<n; k++) { if (fabs(m[j][k]) > largest) { max_row = j; max_column = k; largest = fabs(m[j][k]); } } } if (largest == 0) { /* Found no element larger than 0: the minor is null, */ /* so the matrix is singular. */ free(swap); for(j=0;j<n;j++) free(m[j]); free(m); return(0); } /* Now move the pivot to the upper left corner of the minor. */ if (max_row != i) { for (j=0; j<n; j++) swap[j] = m[i][j]; for (j=0; j<n; j++) m[i][j] = m[max_row][j]; for (j=0; j<n; j++) m[max_row][j] = swap[j]; sign *= -1; } if (max_column != i) { for (j=0; j<n; j++) swap[j] = m[j][i]; for (j=0; j<n; j++) m[j][i] = m[j][max_column]; for (j=0; j<n; j++) m[j][max_column] = swap[j]; sign *= -1; } } /* End of pivoting choice */ /* Now adjust rows i+1...n so as to zero the i-th element */ for (j=i+1; j<n; j++) { if (m[j][i] == 0) continue; /* m[i][i] is the pivot element */ mul_factor = m[j][i]/m[i][i]; /* We cannot trust the subtraction of two very close numbers to return exactly zero because of the roundoff error... */ m[j][i] = 0; for (k=i+1; k<n; k++) { m[j][k] = m[j][k] - mul_factor * m[i][k]; } } } /* Now the determinant is just the products of the diagonal */ /* elements. */ /* We could do this, but a log sum can protect us from */ /* overflows */ /* * result = 1; * * for (i=0;i<n;i++) * { * result *= m[i][i]; * } */ result = 0; for (i=0;i<n;i++) { if (m[i][i] < 0) sign *= -1; result += log(fabs(m[i][i])); } free(swap); for(j=0;j<n;j++) free(m[j]); free(m); return(sign*exp(result)); } int print_matrix(double **matrix, int size) { int i,j; for (i=0;i<size;i++) { for (j=0;j<size;j++) { printf("[%10.5g]",matrix[i][j]); } printf("\n"); } printf("\n"); return(0); }