aaz_plu_r.c

一些矩阵运算的C代码实现：矩阵一般运算

/* gjid03.c      freeware     http://groups.yahoo.com/group/mathc/
*/
/* --------------------------------- INCLUDES ------------------------------- */
#include "y_ahfile.h"
/* ------------------------------------- FUNCTION --------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* -------------------------------------------------------------------------- */
void gjA(int n)
{
int      x;

double pbA[MXR][MXC*C2];mZ mA={MXR,MXC*C2,&pbA[0][0]};
double pbP[MXR][MXC*C2];mZ mP={MXR,MXC*C2,&pbP[0][0]};
double pbL[MXR][MXC*C2];mZ mL={MXR,MXC*C2,&pbL[0][0]};
double pbU[MXR][MXC*C2];mZ mU={MXR,MXC*C2,&pbU[0][0]};
double pbT[MXR][MXC*C2];mZ mT={MXR,MXC*C2,&pbT[0][0]};

double pinvP   [MXR][MXR*C2]; mZ minvP   ={MXR,MXR*C2,&pinvP   [0][0]};
double pinvP_LU[MXR][MXR*C2]; mZ minvP_LU={MXR,MXR*C2,&pinvP_LU[0][0]};

                     mA.rows=n;      mA.cols=n*C2;
                     mP.rows=n;      mP.cols=n*C2;
                     mL.rows=n;      mL.cols=n*C2;
                     mU.rows=n;      mU.cols=n*C2;
                     mT.rows=n;      mT.cols=n*C2;

                  minvP.rows=n;   minvP.cols=n*C2;
               minvP_LU.rows=n;minvP_LU.cols=n*C2;
/*-------------------------------- PROGRAM ---------------------------------- */
  rnd0_mZ(&mA,3);

  clrscrn();
  printf("                                                       A :");
  p_mZ(&mA);
  PLU_mZ(&mA,&mP,&mL,&mU);
  printf("                                                       P :");
  p_mZ(&mP);
  printf("                                                       L :");
  p_mZ(&mL);
  printf("                                                       U :");
  p_mZ(&mU);
  getchar();

  clrscrn();
  printf("                                                 PA = LU :\n\n");
  printf("                                                      LU :");
  mul_mZ(&mL,&mU,&mT);
  p_mZ(&mT);
  printf("                                                      PA :");
  mul_mZ(&mP,&mA,&mT);
  p_mZ(&mT);
  getchar();

  clrscrn();
  printf("                                             A = invP_LU :\n\n");
  printf("                                                       A :");
  p_mZ(&mA);
  printf("                                                 invP_LU :");
  inverse_mZ(&mP,&minvP);
  mul_mZ(&mL,&mU,&mT);
  mul_mZ(&minvP, &mT,&minvP_LU);
  p_mZ(&minvP_LU);
}

void hello(void)
{
  printf("\n Freeware by : xhunga bernard 06000 Nice France.\n");
  printf("\n                                                  ");
  printf("\n           PLU decomposition.                     ");
  printf("\n           *****************                    \n");
  printf("\n The size of mA : 2x2 ... 4x4.                  \n");
  printf("\n           ex : \"LUFR.exe\"                    \n");
  printf("\n You can choose the size :                      \n");
  printf("\n           ex : \"LUFR.exe 2\" -> mA 2x2          ");
  printf("\n                \"LUFR.exe 3\" -> mA 3x3          ");
  printf("\n                \"LUFR.exe 4\" -> mA 4x4          ");
  printf("\n\n                                                ");
  
  printf("\n                                                  ");
  printf("\n Press return to continue                         ");
  printf("\n                                                  ");

  getchar();
}

/* --------------------------------- MAIN ----------------------------------- */
int main(int argc, char *argv[])
{
time_t t;
int    n;

/*-------------------------------- PROGRAM ---------------------------------- */
  srand(time(&t));
  hello();
do
{
if( argc == 2 && atoi(argv[1]))
 {  n = atoi(argv[1]); if( n > 4){n=4;}}
else
 {  n = rndp_I(4);}
if( n < 2){n=2;}


  gjA(4);


printf("\n Press 'X' Return to stop\n");
} while(toupper(getchar())!='X');

  return 0;
}