-
Notifications
You must be signed in to change notification settings - Fork 0
/
CGSerial.cpp
126 lines (110 loc) · 3.79 KB
/
CGSerial.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
/*
* File: CGSerial.cpp
* Author: zhakov
*
* Created on 20 Декабрь 2012 г., 23:48
*/
#include "CGSerial.h"
#include "stdlib.h"
#include "math.h"
#include "stdio.h"
#include "matrixHelpers.h"
CGSerial::CGSerial() {
}
CGSerial::CGSerial(const CGSerial& orig) {
}
CGSerial::~CGSerial() {
}
void CGSerial::resultCalculation(double** pMatrix, double* pVector, double* pResult, int Size) {
double *CurrentApproximation, *PreviousApproximation;
double *CurrentGradient, *PreviousGradient;
double *CurrentDirection, *PreviousDirection;
double *Denom, *tempPointer;
double Step;
int Iter = 1, MaxIter = Size + 1;
float Accuracy = 0.00001f;
//Allocating memory
CurrentApproximation = new double[Size];
PreviousApproximation = new double[Size];
CurrentGradient = new double[Size];
PreviousGradient = new double[Size];
CurrentDirection = new double[Size];
PreviousDirection = new double[Size];
Denom = new double[Size];
for (int i = 0; i < Size; i++) {
PreviousApproximation[i] = 0;
PreviousDirection[i] = 0;
PreviousGradient[i] = -pVector[i];
}
do {
if (Iter > 1) {
tempPointer = PreviousApproximation;
PreviousApproximation = CurrentApproximation;
CurrentApproximation = tempPointer;
//SwapPointers(PreviousApproximation, CurrentApproximation);
tempPointer = PreviousGradient;
PreviousGradient = CurrentGradient;
CurrentGradient = tempPointer;
//SwapPointers(PreviousGradient, CurrentGradient);
tempPointer = PreviousDirection;
PreviousDirection = CurrentDirection;
CurrentDirection = tempPointer;
//SwapPointers(PreviousDirection, CurrentDirection);
}
//compute gradient
//#pragma omp parallel for
for (int i = 0; i < Size; i++) {
CurrentGradient[i] = -pVector[i];
for (int j = 0; j < Size; j++)
CurrentGradient[i] += pMatrix[i][j] * PreviousApproximation[j];
}
//compute direction
double IP1 = 0, IP2 = 0;
//#pragma omp parallel for reduction(+:IP1,IP2)
for (int i = 0; i < Size; i++) {
IP1 += CurrentGradient[i] * CurrentGradient[i];
IP2 += PreviousGradient[i] * PreviousGradient[i];
}
//#pragma omp parallel for
for (int i = 0; i < Size; i++) {
CurrentDirection[i] = -CurrentGradient[i] +
PreviousDirection[i] * IP1 / IP2;
}
//compute size step
IP1 = 0;
IP2 = 0;
//#pragma omp parallel for reduction(+:IP1,IP2)
for (int i = 0; i < Size; i++) {
Denom[i] = 0;
for (int j = 0; j < Size; j++)
Denom[i] += pMatrix[i][j] * CurrentDirection[j];
IP1 += CurrentDirection[i] * CurrentGradient[i];
IP2 += CurrentDirection[i] * Denom[i];
}
Step = -IP1 / IP2;
//#pragma omp parallel for
for (int i = 0; i < Size; i++) {
CurrentApproximation[i] = PreviousApproximation[i] + Step * CurrentDirection[i];
}
Iter++;
} while
((diff(PreviousApproximation, CurrentApproximation, Size) > Accuracy)
&& (Iter < MaxIter));
for (int i = 0; i < Size; i++)
pResult[i] = CurrentApproximation[i];
iterationsCount = Iter;
delete CurrentApproximation;
delete PreviousApproximation;
delete CurrentGradient;
delete PreviousGradient;
delete CurrentDirection;
delete PreviousDirection;
delete Denom;
}
double CGSerial::diff(double *vector1, double* vector2, int Size) {
double sum = 0;
for (int i = 0; i < Size; i++) {
sum += fabs(vector1[i] - vector2[i]);
}
return sum;
}