Obtaining K^-1 (invmatrix) for matrix K , using C++, Visual studio

• Question

• Listing for C99

testinv.cpp

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include "invmatrix.h"

int main (void)
{
int order=4;

double **K =(double **) malloc(order*sizeof(double *));
int i;
for (i = 0; i < order; i++) {
K[i]=(double *) malloc(order * sizeof(double ));
}
double **IK =(double **) malloc(order*sizeof(double *));

for (i = 0; i < order; i++) {
IK[i]=(double *) malloc(order * sizeof(double ));
}

K[0][0]=1;
K[0][1]=1;
K[0][2]=1;
K[0][3]=1;
K[1][0]=1;
K[1][1]=0.6975;
K[1][2]=0.48651;
K[1][3]=0.33934;
K[2][0]=1;
K[2][1]=0.2775;
K[2][2]=0.07701;
K[2][3]=0.02137;
K[3][0]=1;
K[3][1]=0.0002;
K[3][2]=3.9996e-8;
K[3][3]=7.9988e-12;

int m1,n1;

for (m1=0;m1<4;m1++){ //row
for (n1=0;n1<4;n1++){ //col

printf ("  %e ",K[m1][n1] );

}

printf ("\n" );
}

printf ("\n Obtaining K^-1\n");
getInvm(  K , IK ,4  );
for (m1=0; m1<4; m1++){ //row
for (n1=0;n1<4; n1++){ //col
//IK[m1][n1]=n1*m1;
printf ("  %le ", IK[m1][n1] );

}

printf ("\n" );
}

getch();
return 0;
}

invmatrix.h

/*
r  c                         1n
k[m][n]={ 1     1      1       1           }
1   0.6975  0.48651  0.33934
1  0.2775   0.07701  0.02137
m1 1  0.0002   3.9996e-8 7.9988e-12 mn

deta=0.0177

AS[m][n]= -3.144*10^{-6} ; 1.1118*10^{-5}; -2.942*10^{-5};0.0178
0.0157          -0.0556          0.1472         -0.1073
-0.0792         0.2559           -0.358         0.1813
0.0812          -0.2003         0.2109         -0.0918

K^-1=  -0.0002;0.0006;-0.0017;1.0012
0.8868;-3.1355;8.2935;-6.0448
-4.4631;14.4227;-20.1761;10.2164
4.5765;-11.2879;11.8842;-5.1728

*/

double  determinant(long double **matrix, int order);
double algdop(long double **matrix, int order, int i,int j);

double** submatrix( double **matrix, int order, int i, int j)
{

double **subm;
int p, q;    // Indexes for matrix
int a = 0, b;  // Indexes for subm
//subm = new  double* [order - 1];

subm = (double **) malloc((order-1)*sizeof(double *));

for(p = 0; p < order; p++) {
if(p==i) continue;    //Skip ith row

//subm[a] = new  double[order - 1];
subm[a]=(double *) malloc((order-1)*sizeof(double));
b = 0;

for(q = 0; q< order; q++) {
if(q==j) continue;  //Skip jth column
subm[a][b++] = matrix[p][q];
}
a++; //Increment row index
}
return subm;
}

double  determinant( double **matrix, int order)
{
double determ = 0;
if(order == 1)
return **matrix; //Return the element if the matrix is of order one
if(order == 2){
determ=(matrix[0][0]*matrix[1][1]-matrix[0][1]*matrix[1][0]);
return determ;
}
int i;

for(i = 0; i < order; i++)
determ +=  (pow(-1.0,i)) * matrix[i][0] * determinant(submatrix(matrix,order, i, 0), order - 1);
//delete  matrix;
return determ;
}

double algdop( double **matrix, int order, int i,int j)
{
return (pow(-1,i+j))*matrix[i][0] * determinant(submatrix(matrix,order, i, j), order - 1);
}

void getInvm( double **K , double **INVM, int order  )
{

int m,n;
double AS[order][order];
double det;
det=determinant(K, order);
for (m=0;m<order;m++){   // rows
for (n=0;n<order;n++){   // columns
// for a(m,n) A(m,n)-algdop m,n
INVM[m][n] =algdop(K, order, n,m)/det;  //A*={A(n,m)*},INVK=A*/determ(K)

}
}

return ;
}

Saturday, September 26, 2015 7:50 PM

All replies

• Use long double ,le , lu , , cstdio, cstdlib , not use getch () if run  in Visual Studio (C++).
Saturday, September 26, 2015 7:51 PM
• is it good for long double ? (without "new",using malloc for 2D pointered arrays)
Saturday, September 26, 2015 7:53 PM
• Example of solving equation with invmatrix K

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include "invmatrix.h"

int main (void)
{
int order=4;

double **K =(double **) malloc(order*sizeof(double *));
int i;
for (i = 0; i < order; i++) {
K[i]=(double *) malloc(order * sizeof(double ));
}
double **IK =(double **) malloc(order*sizeof(double *));

for (i = 0; i < order; i++) {
IK[i]=(double *) malloc(order * sizeof(double ));
}

K[0][0]=1;
K[0][1]=1;
K[0][2]=1;
K[0][3]=1;
K[1][0]=1;
K[1][1]=0.6975;
K[1][2]=0.48651;
K[1][3]=0.33934;
K[2][0]=1;
K[2][1]=0.2775;
K[2][2]=0.07701;
K[2][3]=0.02137;
K[3][0]=1;
K[3][1]=0.0002;
K[3][2]=3.9996e-8;
K[3][3]=7.9988e-12;

int m1,n1;

for (m1=0;m1<4;m1++){ //row
for (n1=0;n1<4;n1++){ //col

printf ("  %e ",K[m1][n1] );

}

printf ("\n" );
}

printf ("\n Obtaining K^-1\n");
getInvm(  K , IK ,order  );
for (m1=0; m1<4; m1++){ //row
for (n1=0;n1<4; n1++){ //col
//IK[m1][n1]=n1*m1;
printf ("  %le ", IK[m1][n1] );

}

printf ("\n" );
}

double A[order];
double F[order];

F[0]=1.0;
F[1]=0.7015;
F[2]=0.39743;
F[3]=0.25046;

printf ("\n F Matrix \n" );

for (n1=0; n1<order; n1++){ //row

printf (" F[%d]= %le ",n1, F[n1] );

printf ("\n" );
}

for (m1=0; m1<order; m1++)   { //row
A[m1]=0;

for (n1=0;n1<order; n1++){ //col
//A=K^-1*F, A, matrix-column, F-matrix-column

A[m1]=A[m1]+IK[m1][n1]*F[n1];// matrix-column of result

}

}

printf ("\n A Matrix \n" );

for (n1=0; n1<order; n1++){ //row

printf (" A[%d]= %le ",n1, A[n1] );

printf ("\n" );
}

getch();
return 0;
}

Saturday, September 26, 2015 8:45 PM