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++)
for (j=0; j< M3. getactualsize (); j++)
M3. getvalue (i,j,rv,xyz);
// выделение памяти для матрицы результатов
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;
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 ();
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) {
left. getvalue (i,j,leftvalue,success);
right. getvalue (i,j,rightvalue,success);
sum = leftvalue + rightvalue;
void setactualsize (int newactualsize) {
if (newactualsize > maxsize)
maxsize = newactualsize;
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;
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 ();
D value;
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++) {
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++) {
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] ;
sum += data [k*maxsize+j] * ( (i==k)? 1.0: data [i*maxsize+k]);
data [i*maxsize+j] = - sum;
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