Sunday, January 11, 2009

Inverse ,determinant ,adjoint for NxN matrix

Source URL


With slight modification. Number of
rows will be input by user and adjoint
will also be printed out.


#include<stdlib.h>
#include<stdio.h>
#include<assert.h>
#include<errno.h>
#include <iostream.h>

float determinant(float **matrix,
int size,
int * used_rows,
int * used_cols,
int depth);
float coefficient(float **matrix,int size, int row, int col);
void print_matrix(FILE * fptr,float ** mat,int rows, int cols);
float test_data[8][8] = {
{4.0, 2.0, 4.0, 5.0, 4.0, -2.0, 4.0, 5.0},
{4.0, 2.0, 4.0, 5.0, 3.0, 9.0, 12.0, 1.0 },
{3.0, 9.0, -13.0, 15.0, 3.0, 9.0, 12.0, 15.0},
{3.0, 9.0, 12.0, 15.0, 4.0, 2.0, 7.0, 5.0 },
{2.0, 4.0, -11.0, 10.0, 2.0, 4.0, 11.0, 10.0 },
{2.0, 4.0, 11.0, 10.0, 3.0, -5.0, 12.0, 15.0 },
{1.0, -2.0, 4.0, 10.0, 3.0, 9.0, -12.0, 15.0 } ,
{1.0, 2.0, 4.0, 10.0, 2.0, -4.0, -11.0, 10.0 } ,
};
#define ROWS 8
void print_adjoint(FILE* fptr,float ** mat,int rows,int cols,float det)
{
int i,j;
for(i=0;i<rows;i++)
{
for(j=0;j<cols;j++)
{
fprintf(fptr,"%10.4f ",mat[i][j]*det);
}
fprintf(fptr,"\n");
}
fflush(fptr);

}

int main(int argc, char **argv)
{

float **matrix;
float **inverse;
int rows,i,j;
float determ;
int * used_rows, * used_cols;

//rows = ROWS;
cout << "Input number of rows" << endl;
cin >> rows;
cout << "Now enter the input values" << endl;
/* Allocate markers to record rows and columns to be skipped */

/* during determinant calculation */
used_rows = (int *) malloc(rows*sizeof(*used_rows));
used_cols = (int *) malloc(rows*sizeof(*used_cols));

/* Allocate working copy of matrix and initialize it from static copy */
matrix = (float **) malloc(rows*sizeof(*matrix));
inverse = (float **) malloc(rows*sizeof(*inverse));
for(i=0;i<rows;i++)
{
matrix[i] = (float *) malloc(rows*sizeof(**matrix));
inverse[i] = (float *) malloc(rows*sizeof(**inverse));
for(j=0;j<rows;j++)
//matrix[i][j] = test_data[i][j];
cin >> matrix[i][j];
}

/* Compute and print determinant */
printf("The determinant of\n\n");
print_matrix(stdout,matrix,rows,rows);
determ=determinant(matrix,rows,used_rows,used_cols,0);
printf("\nis %f\n",determ);
fflush(stdout);
assert(determ!=0);

for(i=0;i<rows;i++)
{
for(j=0;j<rows;j++)
{
inverse[j][i] = coefficient(matrix,rows,i,j)/determ;
}
}

printf("The inverse is\n\n");
print_matrix(stdout,inverse,rows,rows);
printf("\n\nThe adjoint is\n\n");
print_adjoint(stdout,inverse,rows,rows,determ);

return 0;
}

float determinant(float **matrix,int size, int * used_rows, int * used_cols, int depth)
{
int col1, col2, row1, row2;
int j,k;
float total=0;
int sign = 1;

/* Find the first unused row */
for(row1=0;row1<size;row1++)
{
for(k=0;k<depth;k++)
{
if(row1==used_rows[k]) break;
}
if(k>=depth) /* this row is not used */
break;
}
assert(row1<size);

if(depth==(size-2))
{
/* There are only 2 unused rows/columns left */

/* Find the second unused row */
for(row2=row1+1;row2<size;row2++)
{
for(k=0;k<depth;k++)
{
if(row2==used_rows[k]) break;
}
if(k>=depth) /* this row is not used */
break;
}
assert(row2<size);

/* Find the first unused column */
for(col1=0;col1<size;col1++)
{
for(k=0;k<depth;k++)
{
if(col1==used_cols[k]) break;
}
if(k>=depth) /* this column is not used */
break;
}
assert(col1<size);

/* Find the second unused column */
for(col2=col1+1;col2<size;col2++)
{
for(k=0;k<depth;k++)
{
if(col2==used_cols[k]) break;
}
if(k>=depth) /* this column is not used */
break;
}
assert(col2<size);

/* Determinant = m11*m22-m12*m21 */
return matrix[row1][col1]*matrix[row2][col2]-matrix[row2][col1]*matrix[row1][col2];
}

/* There are more than 2 rows/columns in the matrix being processed */
/* Compute the determinant as the sum of the product of each element */
/* in the first row and the determinant of the matrix with its row */
/* and column removed */
total = 0;

used_rows[depth] = row1;
for(col1=0;col1<size;col1++)
{
for(k=0;k<depth;k++)
{
if(col1==used_cols[k]) break;
}
if(k<depth) /* This column is used - skip it */
continue;
used_cols[depth] = col1;
total += sign*matrix[row1][col1]*determinant(matrix,size,used_rows,used_cols,depth+1);
sign=(sign==1)?-1:1;
}
return total;

}

void print_matrix(FILE * fptr,float ** mat,int rows, int cols)
{
int i,j;
for(i=0;i<rows;i++)
{
for(j=0;j<cols;j++)
{
fprintf(fptr,"%10.4f ",mat[i][j]);
}
fprintf(fptr,"\n");
}
fflush(fptr);
}


float coefficient(float **matrix,int size, int row, int col)
{
float coef;
int * ur, *uc;

ur = (int*)malloc(size*sizeof(matrix));
uc = (int*)malloc(size*sizeof(matrix));
ur[0]=row;
uc[0]=col;
coef = (((row+col)%2)?-1:1)*determinant(matrix,size,ur,uc,1);
return coef;
}

No comments:

Blog Archive