Смекни!
smekni.com

Алгоритм компактного хранения и решения СЛАУ высокого порядка (стр. 4 из 5)

          // Create 8*Tetraedr(4)

          Convert1024(FE,NumFE);

        }

      Convert(X,Y,Z,FE,NumFE,Bounds,NumBounds);

       return!Output(fn3,X,Y,Z,FE,NumPoints,0,NumFE,Bounds,NumBounds);

    }

   if (CurrentType == BASE3D_8)

    {

       if (!ReadCubeData(fn1,fn2,in1,in2,X,Y,Z,FE,NumPoints,NumFE))return false;

       if (!strcmp(Op,"/6"))

        {

          // Create 6*Tetraedr(4)

          Convert824(FE,NumFE);

        }

      Convert(X,Y,Z,FE,NumFE,Bounds,NumBounds);

       return!Output(fn3,X,Y,Z,FE,NumPoints,0,NumFE,Bounds,NumBounds);

    }

  return false;

}

void Convert824(Matrix<DWORD>&FE,DWORD& NumFE)

{

  Matrix<DWORD> nFE(6 * NumFE,4);

  DWORD  data[][4] = {

                      { 0,2,3,6 },

                      { 4,5,1,7 },

                      { 0,4,1,3 },

                      { 6,7,3,4 },

                      { 1,3,7,4 },

                      { 0,4,6,3 }

                     },

         Current = 0;

  for (DWORD i = 0; i < NumFE; i++)

   for (DWORD j = 0; j < 6; j++)

    {

      for (DWORD k = 0; k < 4; k++)

      nFE[Current][k] = FE[i][data[j][k]];

      Current++;

    }

  CurrentType = BASE3D_4;

  NumFE = Current;

  FE    = nFE;

}

void Convert1024(Matrix<DWORD>&FE,DWORD& NumFE)

{

  Matrix<DWORD> nFE(8 * NumFE,4);

  DWORD  data[][4] = {

                      { 3,7,8,9 },

                      { 0,4,7,6 },

                      { 2,5,9,6 },

                      { 7,9,8,6 },

                      { 4,8,5,1 },

                      { 4,5,8,6 },

                      { 7,8,4,6 },

                      { 8,9,5,6 }

                     },

         Current = 0;

  for (DWORD i = 0; i < NumFE; i++)

   for (DWORD j = 0; j < 8; j++)

    {

      for (DWORD k = 0; k < 4; k++)

      nFE[Current][k] = FE[i][data[j][k]];

      Current++;

    }

  CurrentType = BASE3D_4;

  NumFE = Current;

  FE    = nFE;

}

bool ReadTetraedrData(char* fn1,char*fn2,FILE* in1,FILE* in2,Vector<double>& X,Vector<double>&Y,Vector<double>& Z,

                     Matrix<DWORD>& FE,DWORD& NumPoints,DWORD& NumFE)

{

   double         tx,

                  ty,

                  tz;

   char           TextBuffer[1001];

   DWORD          Current = 0L,

                  tmp;

   while (!feof(in1))

    {

      if (fgets(TextBuffer,1000,in1) ==NULL)

       {

         if (feof(in1)) break;

         printf("Unable read file%s",fn1);

         fclose(in1);

         fclose(in2);

         return false;

       }

      if (sscanf(TextBuffer,"%ld %lf%lf %lf", &NumPoints,&tx,&ty,&tz) != 4) continue;

      X[Current] = tx;

      Y[Current] = ty;

      Z[Current] = tz;

      if (++Current == 999)

       {

         Vector<double> t1 = X,

                        t2 = Y,

                        t3 = Z;

         X.ReSize(2 * X.Size());

         Y.ReSize(2 * X.Size());

         Z.ReSize(2 * X.Size());

         for (DWORD i = 0; i < Current;i++)

          {

            X[i] = t1[i];

            Y[i] = t2[i];

            Z[i] = t3[i];

          }

       }

      if (Current % 100 == 0)

       printf("Line:%ld&bsol;r",Current);

    }

   fclose(in1);

   printf("                           &bsol;r");

   NumPoints = Current;

   Current = 0L;

   while (!feof(in2))

    {

      if (fgets(TextBuffer,1000,in2) ==NULL)

       {

         if (feof(in2)) break;

         printf("Unable read file%s",fn2);

         fclose(in2);

         return false;

       }

      if (sscanf(TextBuffer,"%d %d %d%d %d %ld %ld %ld %ld %ld %ld %ld %ld",

         &tmp,&tmp,&tmp,&tmp,&tmp,

         &FE[Current][0],&FE[Current][1],&FE[Current][2],&FE[Current][3],

         &FE[Current][4],&FE[Current][5],&FE[Current][6],&FE[Current][7])!= 13) continue;

      if (fgets(TextBuffer,1000,in2) ==NULL)

       {

         printf("Unable read file%s",fn2);

         fclose(in2);

         return false;

       }

      if (sscanf(TextBuffer,"%ld%ld",&FE[Current][8],&FE[Current][9]) != 2)

       {

         printf("Unable read file%s",fn2);

         fclose(in2);

         return false;

       }

{

       if (fabs((tx = 0.5*(X[FE[Current][0]- 1] + X[FE[Current][1] - 1])) - X[FE[Current][4] - 1]) > Eps)

        X[FE[Current][4] - 1] = tx;

       if (fabs((ty = 0.5*(Y[FE[Current][0]- 1] + Y[FE[Current][1] - 1])) - Y[FE[Current][4] - 1]) > Eps)

          Y[FE[Current][4] - 1] = ty;

       if (fabs((tz = 0.5*(Z[FE[Current][0]- 1] + Z[FE[Current][1] - 1])) - Z[FE[Current][4] - 1]) > Eps)

          Z[FE[Current][4] - 1] = tz;

       if (fabs((tx = 0.5*(X[FE[Current][2]- 1] + X[FE[Current][1] - 1])) - X[FE[Current][5] - 1]) > Eps)

          X[FE[Current][5] - 1] = tx;

       if (fabs((ty = 0.5*(Y[FE[Current][2]- 1] + Y[FE[Current][1] - 1])) - Y[FE[Current][5] - 1]) > Eps)

          Y[FE[Current][5] - 1] = ty;

       if (fabs((tz = 0.5*(Z[FE[Current][2]- 1] + Z[FE[Current][1] - 1])) - Z[FE[Current][5] - 1]) > Eps)

          Z[FE[Current][5] - 1] = tz;

       if (fabs((tx = 0.5*(X[FE[Current][0]- 1] + X[FE[Current][2] - 1])) - X[FE[Current][6] - 1]) > Eps)

          X[FE[Current][6] - 1] = tx;

       if (fabs((ty = 0.5*(Y[FE[Current][0]- 1] + Y[FE[Current][2] - 1])) - Y[FE[Current][6] - 1]) > Eps)

          Y[FE[Current][6] - 1] = ty;

       if (fabs((tz = 0.5*(Z[FE[Current][0]- 1] + Z[FE[Current][2] - 1])) - Z[FE[Current][6] - 1]) > Eps)

          Z[FE[Current][6] - 1] = tz;

       if (fabs((tx = 0.5*(X[FE[Current][0]- 1] + X[FE[Current][3] - 1])) - X[FE[Current][7] - 1]) > Eps)

          X[FE[Current][7] - 1] = tx;

       if (fabs((ty = 0.5*(Y[FE[Current][0]- 1] + Y[FE[Current][3] - 1])) - Y[FE[Current][7] - 1]) > Eps)

          Y[FE[Current][7] - 1] = ty;

       if (fabs((tz = 0.5*(Z[FE[Current][0]- 1] + Z[FE[Current][3] - 1])) - Z[FE[Current][7] - 1]) > Eps)

          Z[FE[Current][7] - 1] = tz;

       if (fabs((tx = 0.5*(X[FE[Current][3]- 1] + X[FE[Current][1] - 1])) - X[FE[Current][8] - 1]) > Eps)

          X[FE[Current][8] - 1] = tx;

       if (fabs((ty = 0.5*(Y[FE[Current][3]- 1] + Y[FE[Current][1] - 1])) - Y[FE[Current][8] - 1]) > Eps)

          Y[FE[Current][8] - 1] = ty;

       if (fabs((tz = 0.5*(Z[FE[Current][3]- 1] + Z[FE[Current][1] - 1])) - Z[FE[Current][8] - 1]) > Eps)

          Z[FE[Current][8] - 1] = tz;

       if (fabs((tx = 0.5*(X[FE[Current][3]- 1] + X[FE[Current][2] - 1])) - X[FE[Current][9] - 1]) > Eps)

          X[FE[Current][9] - 1] = tx;

       if (fabs((ty = 0.5*(Y[FE[Current][3]- 1] + Y[FE[Current][2] - 1])) - Y[FE[Current][9] - 1]) > Eps)

          Y[FE[Current][9] - 1] = ty;

       if (fabs((tz = 0.5*(Z[FE[Current][3]- 1] + Z[FE[Current][2] - 1])) - Z[FE[Current][9] - 1]) > Eps)

          Z[FE[Current][9] - 1] = tz;

}

      if (++Current == 999)

       {

         Matrix<DWORD> t = FE;

         FE.ReSize(2 * FE.Size1(),10);

         for (DWORD i = 0; i < Current;i++)

          for (DWORD j = 0; j < 10; j++)

           FE[i][j] = t[i][j];

       }

      if (Current % 100 == 0)

       printf("Line:%ld&bsol;r",Current);

    }

   NumFE = Current;

   for (DWORD i = 0; i < NumFE; i++)

    for (DWORD j = 0; j < 10; j++)

      FE[i][j]--;

   printf("                           &bsol;r");

   return true;

}

bool ReadCubeData(char* fn1,char*fn2,FILE*in1,FILE* in2,Vector<double>& X,Vector<double>&Y,Vector<double>& Z,

                     Matrix<DWORD>& FE,DWORD& NumPoints,DWORD& NumFE)

{

   double         tx,

                  ty,

                  tz;

   char           TextBuffer[1001];

   DWORD          Current = 0L,

                  tmp;

   while (!feof(in1))

    {

      if (fgets(TextBuffer,1000,in1) ==NULL)

       {

         if (feof(in1)) break;

         printf("Unable read file%s",fn1);

         fclose(in1);

         fclose(in2);

         return false;

       }

      if (sscanf(TextBuffer,"%ld %lf%lf %lf", &NumPoints,&tx,&ty,&tz) != 4) continue;

      X[Current] = tx;

      Y[Current] = ty;

      Z[Current] = tz;

      if (++Current == 999)

       {

         Vector<double> t1 = X,

                        t2 = Y,

                        t3 = Z;

         X.ReSize(2 * X.Size());

         Y.ReSize(2 * X.Size());

         Z.ReSize(2 * X.Size());

         for (DWORD i = 0; i < Current;i++)

          {

            X[i] = t1[i];

            Y[i] = t2[i];

            Z[i] = t3[i];

          }

       }

      if (Current % 100 == 0)

       printf("Line:%ld&bsol;r",Current);

    }

   fclose(in1);

   printf("                           &bsol;r");

   NumPoints = Current;

   Current = 0L;

   while (!feof(in2))

    {

      if (fgets(TextBuffer,1000,in2) ==NULL)

       {

         if (feof(in2)) break;

         printf("Unable read file%s",fn2);

         fclose(in2);

         return false;

       }

      if (sscanf(TextBuffer,"%d %d %d%d %d %ld %ld %ld %ld %ld %ld %ld %ld",

         &tmp,&tmp,&tmp,&tmp,&tmp,

         &FE[Current][0],&FE[Current][1],&FE[Current][3],&FE[Current][2],

         &FE[Current][4],&FE[Current][5],&FE[Current][7],&FE[Current][6])!= 13) continue;

      if (++Current == 999)

       {

         Matrix<DWORD> t = FE;

         FE.ReSize(2 * FE.Size1(),10);

         for (DWORD i = 0; i < Current;i++)

          for (DWORD j = 0; j < 10; j++)

           FE[i][j] = t[i][j];

       }

      if (Current % 100 == 0)

       printf("Line:%ld&bsol;r",Current);}

   NumFE = Current;

   for (DWORD i = 0; i < NumFE; i++)

    for (DWORD j = 0; j < 10; j++)

      FE[i][j]--;

   printf("                           &bsol;r");

return true;}ПРИЛОЖЕНИЕ 2.

Исходный текст программы, реализующей алгоритмкомпактного хранения и решения СЛАУ высокого порядка.

#include "matrix.h"

class RVector

{

  private:

   Vector<double> Buffer;

  public:

   RVector(void) {}

  ~RVector() {}

   RVector(DWORD Size) {Buffer.ReSize(Size);  }

   RVector(RVector& right) { Buffer =right.Buffer;  }

   RVector(Vector<double>& right){ Buffer = right;  }

   DWORD Size(void)    { returnBuffer.Size(); }

   void  ReSize(DWORD Size) {Buffer.ReSize(Size); }

   double&  operator [] (DWORD i) {return Buffer[i]; }

   RVector& operator = (RVector&right) { Buffer = right.Buffer; return *this; }

   RVector& operator =(Vector<double>& right) { Buffer = right; return *this; }

   void Sub(RVector&);

   void Sub(RVector&,double);

   void Mul(double);

   void Add(RVector&);

   friend doubleNorm(RVector&,RVector&);

};

class TSMatrix

{

  private:

    Vector<double>   Right;

    Vector<double>*  Array;

    Vector<DWORD>*   Links;

    uint             Dim;

    DWORD            Size;

  public:

   TSMatrix(void) { Size = 0; Dim = 0; Array= NULL; Links = NULL; }

  TSMatrix(Vector<DWORD>*,DWORD,uint);

  ~TSMatrix(void) { if (Array) delete []Array; }

   Vector<double>& GetRight(void){ return Right; }

   DWORD GetSize(void) { return Size; }

   uint  GetDim(void) { return Dim; }

   Vector<double>& GetVector(DWORDi) { return Array[i]; }

   Vector<DWORD>*  GetLinks(void) {return Links; }

   void  SetLinks(Vector<DWORD>* l) {Links = l; }

   void Add(Matrix<double>&,Vector<DWORD>&);

   void Add(DWORD I, DWORD L, DWORD J, DWORDK, double v)

   {

     DWORD Row = I,

           Col = L * Links[I].Size() * Dim +Find(I,J) * Dim + K;

     Array[Row][Col] += v;

   }

   void Add(DWORD I, double v)

   {

     Right[I] += v;

   }

   DWORD Find(DWORD,DWORD);

   void  Restore(Matrix<double>&);

   void  Set(DWORD,DWORD,double,bool);

   void  Set(DWORD Index1,DWORDIndex2,double value)

   {

     DWORD I   = Index1 / Dim,

           L   = Index1 % Dim,

           J   = Index2 / Dim,

           K   = Index2 % Dim,

           Pos = Find(I,J),

           Row = I,

           Col;

     if (Pos == DWORD(-1)) return;

     Col = L * Links[I].Size() * Dim +Find(I,J) * Dim + K;

     Array[Row][Col] = value;

   }

   bool  Get(DWORD Index1,DWORDIndex2,double& value)

   {

     DWORD I   = Index1 / Dim,

           L   = Index1 % Dim,

           J   = Index2 / Dim,

           K   = Index2 % Dim,

           Pos = Find(I,J),

           Row = I,

           Col;

     value = 0;

     if (Pos == DWORD(-1)) return false;

     Col = L * Links[I].Size() * Dim +Find(I,J) * Dim + K;

     value = Array[Row][Col];

     return true;

   }

   void  Mul(RVector&,RVector&);

   double  Mul(DWORD,RVector&);

   void  write(ofstream&);

   void  read(ifstream&);

};

class RMatrix

{

  private:

    Vector<double> Buffer;