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

Страница 7

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<double> 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<DWORD>* links, DWORD size, uint dim)

{

Dim = dim;

Links = links;

Size = size;

Right.ReSize(Dim * Size);

Array = new Vector<double>[Size];

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

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

}

void TSMatrix::Add(Matrix<double>& FEMatr,Vector<DWORD>& 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<double>& 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<double>[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));

}