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

• Question

• 

  

  0

  Sign in to vote

  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