Розробка алгоритму та програми чисельного розв'язку систем лінійних алгебраїчних рівнянь з розрідженою матрицею

1. Зенкевич О., Морган К. Конечные методы и аппроксимация. - М.: Мир, 1980. - 479 с.

2. Зенкевич О. Метод конечных элементов. - М.: Мир, 1975. - 217 с.

3. Стрэнг Г., Фикс Дж. Теория метода конечных элементов. - М.: Мир, 1977. - 346 с.

4. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: Наука, 1987. - 632 с.

5. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. - М.: Наука, 1984. - 176 с.

6. Бахвалов Н.С. Численные методы. - М.: Наука, 1975. - 153 с.

7. Годунов С.К. Решение систем линейных уравнений. - Новосибирск: Наука, 1980. - 81 с.

8. Gustavson F.G. Some basic techniques for solving sparse matrix algorithms. - New York: Plenum Press, 1972. - 215 p.

9. Джордж А., Лиу Дж. Численное решение больших разрежённых систем уравнений. - М.: Мир, 1984. - 333 с.

10. Rose D.J. A graph theoretic study of the numerical solution of sparse positive definite system of linear equations. - New York: Academic Press, 1972. - 473 p.

11. Мосаковский В.И., Гудрамович В.С., Макеев Е.М. Контактные задачи теории оболочек и стержней. - М.: Машиностроение, 1978. - 718 с.

12. Рвачев В.Л., Шевченко А.Н. Проблемно-ориентированные языки и системы для инженерных расчётов. - К.: Техніка, 1988. - 198 с.

13. Голуб Дж., Ван Лоун Ч. Матричные вычисления. - М.: Мир, 1999. - 548 с.

14. Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. - М.: Наука, 1986. - 288 с.

15. Тьюарсон Р. Разрежённые матрицы. - М.: Мир, 1977. - 191 с.

16. Брамеллер А., Аллан Р., Хэмэм Я. Слабозаполненные матрицы. - М.: Энергия, 1979. - 192 с.

17. Эстербю О., Златев 3. Прямые методы для разрежённых матриц. - М.: Мир, 1987. - 120 с.

18. Писсанецки С. Технология разреженных матриц. - М.: Мир, 1988. - 410 с.

19. Сабоннадьер Ж.-К., Кулон Ж.-Л. Метод конечных элементов и САПР. - М.: Мир, 1989. - 192 с.

20. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. - М.: Наука, 1978. - 592 с.

21. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. - М.: Мир, 1991. - 347 с.

22. Валях Е. Последовательно-параллельные вычисления. - М.: Мир, 1985. - 216 с.

23. Таненбаум Э., Ван Стеен М. Распределённые системы. Принципы и парадигмы. - СПб.: Питер, 2003. - 481 с.

Додаток А

Вихідний текст програми розв'язку СЛАР з розрідженою матрицею

#include <math.h>

#include <stdio.h>

#include <string.h>

#include <stdlib.h>

#include <fstream.h>

#include "matrix.h"

#include "smatrix.h"

#define BASE3D_4 4

#define BASE3D_8 8

#define BASE3D_10 10

const double Eps=1.0E-10;

DWORD CurrentType=BASE3D_10;

class RVector

{

private:

Vector<double> Buffer;

public:

RVector() {}

~RVector() {}

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

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

RVector(Vector<double> &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<double> &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<double> Right;

Vector<double> * Array;

Vector<DWORD> * Links;

int Dim;

DWORD Size;

public:

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

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

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

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

DWORD GetSize(void) { return Size; }

int GetDim(void) { return Dim; }

Vector<double> &GetVector(DWORD i) { 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, 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<double>&);

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]; }

};

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, int 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[INT(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);

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

NLink=j%Dim;

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));

}

bool Output(char * fname, Vector<double> &x, Vector<double> &y, Vector<double> &z, Matrix<DWORD> &tr, DWORD n, DWORD NumNewPoints, DWORD ntr, Matrix<DWORD> &Bounds, DWORD CountBn)

{

char * Label="NTRout";

int type=CurrentType, type1=(type==BASE3D_4 || type==BASE3D_10) ? 3 : 4;

DWORD NewSize, i, j;

ofstream Out;

if (type==BASE3D_4) type1=3;

else

if (type==BASE3D_8) type1=4;

else

if (type==BASE3D_10) type1=6;

Out.open(fname, ios::out | ios::binary);

if (Out.bad()) return true;

Out.write((const char *) Label, 6*sizeof(char));

if (Out.fail()) return true;

Out.write((const char *) &type, sizeof(int));

if (Out.fail()) return true;

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

if (Out.fail())

{

Out.close();

return true;

}

Out.write((const char *) &(NewSize=n+NumNewPoints), sizeof(DWORD));

if (Out.fail()) return true;

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

if (Out.fail())

{

Out.close();

return true;

}

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

{

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

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

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

if (Out.fail())

{

Out.close();

return true;

}

}

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

{

Out.write((const char *) &x[n+i], sizeof(double));

Out.write((const char *) &y[n+i], sizeof(double));

if (Out.fail())

{

Out.close();

return true;

}

}

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

if (Out.fail())

{

Out.close();

return true;

}

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

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

{

DWORD out=tr[i][j];

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

if (Out.fail())

{

Out.close();

return true;

}

}

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

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

{

DWORD out=Bounds[i][j];

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

if (Out.fail())

{

Out.close();

return true;

}

}

}

bool Test(DWORD * a, DWORD * b)

{

bool result;

int NumPoints=3;

if (CurrentType==BASE3D_8) NumPoints=4;

else

if (CurrentType==BASE3D_10) NumPoints=6;

for (int i=0; i<NumPoints; i++)

{

result=false;

for (int j=0; j<NumPoints; j++)

if (b[j]==a[i])

{

result=true;

break;

}

if (result==false) return false;

}

return true;

}

void Convert(Vector<double> &X, Vector<double> &Y, Vector<double> &Z, Matrix<DWORD> &FE, DWORD NumTr, Matrix<DWORD> &Bounds, DWORD &BnCount)

{

int cData8[6][5]={{0, 4, 5, 1, 7}, {6, 2, 3, 7, 0}, {4, 6, 7, 5, 0}, {2, 0, 1, 3, 5}, {1, 5, 7, 3, 4}, {6, 4, 0, 2, 1}}, cData4[4][4]={{0, 1, 2, 3}, {1, 3, 2, 0}, {3, 0, 2, 1}, {0, 3, 1, 2}}, cData10[4][7]={{0, 1, 2, 4, 5, 6, 3}, {0, 1, 3, 4, 8, 7, 2}, {1, 3, 2, 8, 9, 5, 0}, {0, 2, 3, 6, 9, 7, 1}}, cData[6][7], Data[6], l, Num1, Num2, m;

DWORD i, j, p[6], pp[6], Index;

Matrix<DWORD> BoundList(4*NumTr, 6);

double cx, cy, cz, x1, y1, z1, x2, y2, z2, x3, y3, z3;

Bounds.ReSize(4*NumTr, 6);

switch (CurrentType)

{

case BASE3D_4:

Num1=4;

Num2=3;

for (l=0; l<Num1; l++)

for (m=0; m<=Num2; m++) cData[l][m]=cData4[l][m];

break;

case BASE3D_8:

Num1=6;

Num2=4;

for (l=0; l<Num1; l++)

for (m=0; m<=Num2; m++) cData[l][m]=cData8[l][m];

break;

case BASE3D_10:

Num1=4;

Num2=6;

for (l=0; l<Num1; l++)

for (m=0; m<=Num2; m++) cData[l][m]=cData10[l][m];

}

printf("Create bounds...\r");

for (i=0; i<NumTr-1; i++)

for (int j=0; j<Num1; j++)

if (!BoundList[i][j])

{

for (l=0; l<Num2; l++) p[l]=FE[i][cData[j][l]];

for (DWORD k=i+1; k<NumTr; k++)

for (int m=0; m<Num1; m++)

if (!BoundList[k][m])

{

for (int l=0; l<Num2; l++)

pp[l]=FE[k][cData[m][l]];

if (Test(p, pp))

BoundList[i][j]=BoundList[k][m]=1;

}

}

x1=X[FE[i][cData[j][0]]];

y1=Y[FE[i][cData[j][0]]];

z1=Z[FE[i][cData[j][0]]];

x2=X[FE[i][cData[j][1]]];

y2=Y[FE[i][cData[j][1]]];

z2=Z[FE[i][cData[j][1]]];

x3=X[FE[i][cData[j][2]]];

y3=Y[FE[i][cData[j][2]]];

z3=Z[FE[i][cData[j][2]]];

for (l=0; l<Num2; l++) Data[l]=cData[j][l];

if (((cx-x1)*(y2-y1)*(z3-z1)+(cy-y1)*(z2-z1)*(x3-x1)+(y3-y1)*(cz-z1)*(x2-x1)-(x3-x1)*(y2-y1)*(cz-z1)-(y3-y1)*(z2-z1)*(cx-x1)-(cy-y1)*(z3-z1)*(x2-x1))>0)

{

if (CurrentType==BASE3D_4)

{

Data[0]=cData[j][0];

Data[1]=cData[j][2];

Data[2]=cData[j][1];

}

else

if (CurrentType==BASE3D_10)

{

Data[0]=cData[j][0];

Data[1]=cData[j][2];

Data[2]=cData[j][1];

Data[3]=cData[j][5];

Data[5]=cData[j][3];

}

else

{

Data[0]=cData[j][0];

Data[1]=cData[j][3];

Data[2]=cData[j][2];

Data[3]=cData[j][1];

}

}

for (l=0; l<Num2; l++) Bounds[BnCount][l]=FE[i][Data[l]];

BnCount++;

}

void main(int argc,char * argv)

{

char * input1, * input2, * input3, * op="", * sw;

if (argc<5 || argc>6) return;

sw=argv[1];

input1=argv[2];

input2=argv[3];

input3=argv[4];

if (!strcmp(sw, "-t10")) CurrentType=BASE3D_10;

else

if (!strcmp(sw, "-c8")) CurrentType=BASE3D_8;

else

{

printf("Unknown switch %s\n\n", sw);

return;

}

if (argc==6)

{

op=argv[5];

if (strcmp(op, "/8") && strcmp(op, "/6"))

{

printf("Unknown options %s\n\n", op);

return;

}

}

if (CreateFile(input1, input2, input3, op)) printf("OK\n");

}

bool CreateFile(char * fn1, char * fn2, char * fn3, char * Op)

{

FILE * in1, * in2, * in3;

Vector<double> X(1000), Y(1000), Z(1000);

DWORD NumPoints, NumFE, NumBounds=0, tmp;

Matrix<DWORD> FE(1000, 10), Bounds;

if ((in1=fopen(fn1, "r"))==NULL)

{

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

return false;

}

if ((in2=fopen(fn2, "r"))==NULL)

{

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

return false;

}

if (CurrentType==BASE3D_10)

{

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

if (!strcmp(Op, "/8")) 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")) 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\r", Current);

}

fclose(in1);

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;

}

}

}

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

}

fclose(in1);

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

}

NumFE=Current;

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

for (DWORD j=0; j<10; j++) FE[i][j]--;

return true;

}