Алгоритм компактного хранения и решения СЛАУ высокого порядка

Vector 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\r",Current);

}

fclose(in1);

printf(" \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 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\r",Current);

}

NumFE = Current;

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

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

FE[i][j]--;

printf(" \r");

return true;

}

bool ReadCubeData(char* fn1,char*fn2,FILE* in1,FILE*

in2,Vector& X,Vector& Y,Vector& Z,

Matrix& 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 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\r",Current);

}

fclose(in1);

printf(" \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 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\r",Current);}

NumFE = Current;

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

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

FE[i][j]--;

printf(" \r");

return true;}

ПРИЛОЖЕНИЕ 2.

Исходный текст программы, реализующей алгоритм компактного хранения и

решения СЛАУ высокого порядка.

#include "matrix.h"

class RVector

{

private:

Vector Buffer;

public:

RVector(void) {}

~RVector() {}

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

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

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

DWORD Size(void) { return Buffer.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& right) { Buffer = right; return

*this; }

void Sub(RVector&);

void Sub(RVector&,double);

void Mul(double);

void Add(RVector&);

friend double Norm(RVector&,RVector&);

};

class TSMatrix

{

private:

Vector Right;

Vector* Array;

Vector* Links;

uint Dim;

DWORD Size;

public:

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

TSMatrix(Vector*,DWORD,uint);

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

Vector& GetRight(void) { return Right; }

DWORD GetSize(void) { return Size; }

uint GetDim(void) { return Dim; }

Vector& GetVector(DWORD i) { return Array[i]; }

Vector* GetLinks(void) { return Links; }

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

void Add(Matrix&,Vector&);

void Add(DWORD I, DWORD L, DWORD J, DWORD K, 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&);

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

void Set(DWORD Index1,DWORD Index2,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,DWORD Index2,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 Buffer;

DWORD size;

public:

RMatrix(DWORD sz) { size = sz; Buffer.ReSize(size*(size + 1)*0.5);

}

~RMatrix() {}

DWORD Size(void) { return size; }

double& Get(DWORD i,DWORD j) { return Buffer[(2*size + 1 - i)*0.5*i

+ j - i]; }

};

//************************

#include "smatrix.h"

double Norm(RVector& Left,RVector& Right)

{

double Ret = 0;

for (DWORD i = 0; i < Left.Size(); i++)

Ret += Left[i] * Right[i];

return Ret;

}

void RVector::Sub(RVector& Right)

{

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

(*this)[i] -= Right[i];

}

void RVector::Add(RVector& Right)

{

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

(*this)[i] += Right[i];

}

void RVector::Mul(double koff)

{

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

(*this)[i] *= koff;

}

void RVector::Sub(RVector& Right,double koff)

{

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

(*this)[i] -= Right[i]*koff;

}

TSMatrix::TSMatrix(Vector* links, DWORD size, uint dim)

{

Dim = dim;

Links = links;

Size = size;

Right.ReSize(Dim * Size);

Array = new Vector[Size];

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

Array[i].ReSize(Links[i].Size() * Dim * Dim);

}

void TSMatrix::Add(Matrix& FEMatr,Vector& FE)

{

double Res;

DWORD RRow;

for (DWORD i = 0L; i < FE.Size(); i++)

for (DWORD l = 0L; l < Dim; l++)

for (DWORD j = 0L; j < FE.Size(); j++)

for (DWORD k = 0L; k < Dim; k++)

{

Res = FEMatr[i * Dim + l][j * Dim + k];

if (Res) Add(FE[i],l,FE[j],k,Res);

}

for (DWORD i = 0L; i < FE.Size(); i++)

for (DWORD l = 0L; l < Dim; l++)

{

RRow = FE[UINT(i % (FE.Size()))] * Dim + l;

Res = FEMatr[i * Dim + l][FEMatr.Size1()];

if (Res) Add(RRow,Res);

}

}

DWORD TSMatrix::Find(DWORD I,DWORD J)

{

DWORD i;

for (i = 0; i < Links[I].Size(); i++)

if (Links[I][i] == J) return i;

return DWORD(-1);

}

void TSMatrix::Restore(Matrix& Matr)

{

DWORD i,

j,

NRow,

NPoint,

NLink,

Pos;

Matr.ReSize(Size * Dim,Size * Dim + 1);

for (i = 0; i < Size; i++)

for (j = 0; j < Array[i].Size(); j++)

{

NRow = j / (Array[i].Size() / Dim); // Number

of row

NPoint = (j - NRow * (Array[i].Size() / Dim)) / Dim; // Number

of points

NLink = j % Dim; // Number

of link

Pos = Links[i][NPoint];

Matr[i * Dim + NRow][Pos * Dim + NLink] = Array[i][j];

}

for (i = 0; i < Right.Size(); i++) Matr[i][Matr.Size1()] = Right[i];

}

void TSMatrix::Set(DWORD Index,DWORD Position,double Value,bool Case)

{

DWORD Row = Index,

Col = Position * Links[Index].Size() * Dim +

Find(Index,Index) * Dim + Position,

i;

double koff = Array[Row][Col],

val;

if (!Case)

Right[Dim * Index + Position] = Value;

else

{

Right[Index * Dim + Position] = Value * koff;

for (i = 0L; i < Size * Dim; i++)

if (i != Index * Dim + Position)

{

Set(Index * Dim + Position,i,0);

Set(i,Index * Dim + Position,0);

if (Get(i,Index * Dim + Position,val))

Right[i] -= val * Value;

}

}

}

void TSMatrix::Mul(RVector& Arr,RVector& Res)

{

DWORD i,

j,

NRow,

NPoint,

NLink,

Pos;

Res.ReSize(Arr.Size());

for (i = 0; i < Size; i++)

for (j = 0; j < Array[i].Size(); j++)

{

NRow = j / (Array[i].Size() / Dim);

NPoint = (j - NRow * (Array[i].Size() / Dim)) / Dim;

NLink = j % Dim;

Pos = Links[i][NPoint];

Res[i * Dim + NRow] += Arr[Pos * Dim + NLink] * Array[i][j];

}

}

double TSMatrix::Mul(DWORD Index,RVector& Arr)

{

DWORD j,

I = Index / Dim,

L = Index % Dim,

Start = L * (Array[I].Size() / Dim),

Stop = Start + (Array[I].Size() / Dim),

NRow,

NPoint,

NLink,

Pos;

double Res = 0;

for (j = Start; j < Stop; j++)

{

NRow = j / (Array[I].Size() / Dim);

NPoint = (j - NRow * (Array[I].Size() / Dim)) / Dim;

NLink = j % Dim;

Pos = Links[I][NPoint];

Res += Arr[Pos * Dim + NLink] * Array[I][j];

}

return Res;

}

void TSMatrix::write(ofstream& Out)

{

DWORD ColSize;

Out.write((char*)&(Dim),sizeof(DWORD));

Out.write((char*)&(Size),sizeof(DWORD));

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

{

ColSize = Array[i].Size();

Out.write((char*)&(ColSize),sizeof(DWORD));

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

Out.write((char*)&(Array[i][j]),sizeof(double));

}

for (DWORD i = 0; i < Size * Dim; i++)

Out.write((char*)&(Right[i]),sizeof(double));

}

void TSMatrix::read(ifstream& In)

{

DWORD ColSize;

In.read((char*)&(Dim),sizeof(DWORD));

In.read((char*)&(Size),sizeof(DWORD));

if (Array) delete [] Array;

Array = new Vector[Size];

Right.ReSize(Size * Dim);

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

{

In.read((char*)&(ColSize),sizeof(DWORD));

Array[i].ReSize(ColSize);

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

In.read((char*)&(Array[i][j]),sizeof(double));

}

for (DWORD i = 0; i < Size * Dim; i++)

In.read((char*)&(Right[i]),sizeof(double));

}

-----------------------

Рисунок 2.

7

6

5

4

1

7

6

4

3

2

1

7

5

4

1

7

6

5

4

3

6

4

3

2

[pic]

6

3

2

1

7

6

5

2

1

7

6

5

4

3

2

1

7

6

5

4

3

2

1

7

6

5

4

3

2

1

[pic]

Рис. 5.

1

4

3

2

5

7

6

Рисунок 1.

Страницы: 1, 2, 3