for (j=0; j< Pow. getactualsize (); j++)

{

Pow. setvalue (i,j,0);

}

}

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

{

Pow. setvalue (i,0,0);

Pow. setvalue (N+i,0,0);

Pow. setvalue (2*N+i,0,q/E*C (10, i,0));

}

printf ("\n\n\nM1: ");

for (i=0; i < M1. getactualsize (); i++)

{

printf ("\n%d: ", i);

for (j=0; j< M1. getactualsize (); j++)

{

M1. getvalue (i,j,rv,xyz);

printf ("%.10f ",rv);

}

}

// выделение памяти под матрицы

matrix <double> M2 (3*N,3*N);

matrix <double> M3 (3*N,3*N);

// копирование матрицы коэффициентов

M2. copymatrix (M1);

// обращение матрицы

M1. invert ();

// произведение матриц

M3. settoproduct (M1,M2);

// сравнение полученной единичной матрицы с эталоном единичной матрицы

M3.comparetoidentity ();

printf ("\n\n\ninverse: ");

for (i=0; i < M3. getactualsize (); i++)

{

printf ("\n%d: ", i);

for (j=0; j< M3. getactualsize (); j++)

{

M3. getvalue (i,j,rv,xyz);

printf ("%.10f ",rv);

}

}

for (i=0; i < M1. getactualsize (); i++)

{

printf ("\n%d: ", i);

for (j=0; j< M1. getactualsize (); j++)

{

M1. getvalue (i,j,rv,xyz);

printf ("%.10f ",rv);

}

}

// выделение памяти для матрицы результатов

matrix <double> Ans (3*N,3*N);

Ans. settoproduct (M1,Pow);

printf ("\nanswer \n");

for (i=0; i < Ans. getactualsize (); i++)

{

printf ("%d: ", i);

Ans. getvalue (i,0,rv,xyz);

printf ("%.10f ", rv);

printf ("\n");

}

// вывод результата

W = 0;

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

{

Ans. getvalue (2*N+i,0,rv,xyz);

printf ("!!%.10f%.10f%.10f\n", rv, X3 [i], Y3 [i]);

W += rv * sin_0 (X3 [i],a/2) * sin_0 (Y3 [i],a/2);

}

printf ("W: =%.10f", W);

// вывод времени вычисления программы

endwtime = MPI_Wtime ();

printf ("\nwall clock time =%f\n", endwtime-startwtime);

fflush (stdout);

}

MPI_Finalize ();

return 0;

}

Matrix. h:

#ifndef __mjdmatrix_h

#define __mjdmatrix_h

template <class D> class matrix{

int maxsize;

int actualsize;

D* data;

void allocate () {

delete [] data;

data = new D [maxsize*maxsize] ;

};

matrix () {};

matrix (int newmaxsize) {matrix (newmaxsize,newmaxsize); };

public:

matrix (int newmaxsize, int newactualsize) {

if (newmaxsize <= 0) newmaxsize = 5;

maxsize = newmaxsize;

if ( (newactualsize <= newmaxsize) && (newactualsize>0))

actualsize = newactualsize;

else

actualsize = newmaxsize;

data = 0;

allocate ();

};

~matrix () { delete [] data; };

void settoproduct (matrix& left, matrix& right) {

actualsize = left. getactualsize ();

if (maxsize < left. getactualsize ()) {

maxsize = left. getactualsize ();

allocate ();

}

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

for (int j = 0; j < actualsize; j++) {

D sum = 0.0;

D leftvalue, rightvalue;

bool success;

for (int c = 0; c < actualsize; c++) {

left. getvalue (i,c,leftvalue,success);

right. getvalue (c,j,rightvalue,success);

sum += leftvalue * rightvalue;

}

setvalue (i,j,sum);

}

}

void combine (matrix& left, matrix& right) {

actualsize = left. getactualsize ();

if (maxsize < left. getactualsize ()) {

maxsize = left. getactualsize ();

allocate ();

}

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

for (int j = 0; j < actualsize; j++) {

D sum = 0.0;

D leftvalue, rightvalue;

bool success;

left. getvalue (i,j,leftvalue,success);

right. getvalue (i,j,rightvalue,success);

sum = leftvalue + rightvalue;

setvalue (i,j,sum);

}

}

void setactualsize (int newactualsize) {

if (newactualsize > maxsize)

{

maxsize = newactualsize;

allocate ();

}

if (newactualsize >= 0) actualsize = newactualsize;

};

int getactualsize () { return actualsize; };

void getvalue (int row, int column, D& returnvalue, bool& success) {

if ( (row>=maxsize) || (column>=maxsize)

|| (row<0) || (column<0))

{ success = false;

return; }

returnvalue = data [row * maxsize + column] ;

success = true;

};

bool setvalue (int row, int column, D newvalue)

data [row * maxsize + column] = newvalue;

void dumpMatrixValues () {

bool xyz;

double rv;

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

{

std:: cout << "i=" << i << ": ";

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

{

M. getvalue (i,j,rv,xyz);

std:: cout << rv << " ";

}

std:: cout << std:: endl;

}

};

void comparetoidentity () {

int worstdiagonal = 0;

D maxunitydeviation = 0.0;

D currentunitydeviation;

for (int i = 0; i < actualsize; i++) {

currentunitydeviation = data [i*maxsize+i] - 1.;

if (currentunitydeviation < 0.0) currentunitydeviation *= - 1.;

if (currentunitydeviation > maxunitydeviation) {

maxunitydeviation = currentunitydeviation;

worstdiagonal = i;

}

}

int worstoffdiagonalrow = 0;

int worstoffdiagonalcolumn = 0;

D maxzerodeviation = 0.0;

D currentzerodeviation;

for (int i = 0; i < actualsize; i++) {

for (int j = 0; j < actualsize; j++) {

if (i == j) continue;

currentzerodeviation = data [i*maxsize+j] ;

if (currentzerodeviation < 0.0) currentzerodeviation *= - 1.0;

if (currentzerodeviation > maxzerodeviation) {

maxzerodeviation = currentzerodeviation;

worstoffdiagonalrow = i;

worstoffdiagonalcolumn = j;

}

}

}

printf ("Worst diagonal value deviation from unity:%0.5f at row/column%0.3f\n", maxunitydeviation, worstdiagonal);

printf ("Worst off-diagonal value deviation from zero:%0.5f at row%0.3f, column%0.3f\n", maxzerodeviation, worstoffdiagonalrow,worstoffdiagonalcolumn);

};

void copymatrix (matrix& source) {

actualsize = source. getactualsize ();

if (maxsize < source. getactualsize ()) {

maxsize = source. getactualsize ();

allocate ();

}

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

for (int j = 0; j < actualsize; j++) {

D value;

bool success;

source. getvalue (i,j,value,success);

data [i*maxsize+j] = value;

}

};

void invert () {

if (actualsize <= 0) return;

if (actualsize == 1) return;

for (int i=1; i < actualsize; i++) data [i] /= data [0] ;

for (int i=1; i < actualsize; i++) {

for (int j=i; j < actualsize; j++) {

D sum = 0.0;

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

sum += data [j*maxsize+k] * data [k*maxsize+i] ;

data [j*maxsize+i] - = sum;

}

if (i == actualsize-1) continue;

for (int j=i+1; j < actualsize; j++) {

D sum = 0.0;

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

sum += data [i*maxsize+k] *data [k*maxsize+j] ;

data [i*maxsize+j] =

(data [i*maxsize+j] -sum) / data [i*maxsize+i] ;

}

}

for (int i = 0; i < actualsize; i++) // invert L

for (int j = i; j < actualsize; j++) {

D x = 1.0;

if (i! = j) {

x = 0.0;

for (int k = i; k < j; k++)

x - = data [j*maxsize+k] *data [k*maxsize+i] ;

}

data [j*maxsize+i] = x / data [j*maxsize+j] ;

}

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

for (int j = i; j < actualsize; j++) {

if (i == j) continue;

D sum = 0.0;

for (int k = i; k < j; k++)

sum += data [k*maxsize+j] * ( (i==k)? 1.0: data [i*maxsize+k]);

data [i*maxsize+j] = - sum;

}

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

for (int j = 0; j < actualsize; j++) {

D sum = 0.0;

for (int k = ( (i>j)? i: j); k < actualsize; k++)

sum += ( (j==k)? 1.0: data [j*maxsize+k]) *data [k*maxsize+i] ;

data [j*maxsize+i] = sum;

}

};

};

#endif

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