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,doublekoff)
{
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);
}
voidTSMatrix::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,DWORDPosition,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(DWORDIndex,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));
}
Список литературы
ЗенкевичО., Морган К. Конечные методы и аппроксимация // М.: Мир, 1980
ЗенкевичО., Метод конечных элементов // М.: Мир., 1975
СтрэнгГ., Фикс Дж. Теория метода конечных элементов // М.: Мир, 1977
БахваловН.С.,Жидков Н.П., Кобельков Г.М. Численные методы // М.: наука, 1987
ВоеводинВ.В., Кузнецов Ю.А. Матрицы и вычисления // М.:Наука, 1984
БахваловН.С. Численные методы // М.: Наука, 1975
ГодуновС.К. Решение систем линейных уравнений // Новосибирск: Наука, 1980
ГоменюкС.И., Толок В.А. Инструментальная система анализа задач механики деформируемоготвердого тела // Приднiпровськийнауковий вiсник – 1997. – №4.
F.G. Gustavson, “Some basic techniques for solving sparse matrixalgorithms”, // editer by D.J. Rose and R.A.Willoughby, Plenum Press, New York,1972
А.Джордж,Дж. Лиу, Численное решение больших разреженных систем уравнений // Москва, Мир,1984
D.J. Rose, “A graph theoretic study of the numerical solution ofsparse positive definite system of linear equations” // New York, AcademicPress, 1972
МосаковскийВ.И., Гудрамович В.С., Макеев Е.М., Контактные задачи теории оболочек истержней // М.:”Машиностроение”, 1978