parallel gauss, parallel CG added, but tests failed

zhakov · zhakov · commit f6263f32d593 · 2012-12-13T22:18:45.000+06:00
diff --git a/CGParallel.cpp b/CGParallel.cpp
@@ -0,0 +1,120 @@
+/* 
+ * File:   CGParallel.cpp
+ * Author: zhakov
+ * Решение СЛАУ методом сопряжённых градиентов
+ * 
+ * Created on 13 Декабрь 2012 г., 21:37
+ */
+
+#include "CGParallel.h"
+#include "stdlib.h"
+#include "math.h"
+
+CGParallel::CGParallel() {
+}
+
+CGParallel::CGParallel(const CGParallel& orig) {
+}
+
+CGParallel::~CGParallel() {
+}
+
+void CGParallel::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.0001f;
+    //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;
+        //compute new approximation
+#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];
+
+
+    delete CurrentApproximation;
+    delete PreviousApproximation;
+    delete CurrentGradient;
+    delete PreviousGradient;
+    delete CurrentDirection;
+    delete PreviousDirection;
+    delete Denom;
+}
+
+double CGParallel::diff(double *vector1, double* vector2, int Size) {
+    double sum = 0;
+    
+    for (int i = 0; i < Size; i++) {
+        sum += abs(vector1[i] - vector2[i]);
+    }
+    
+    return sum;
+}
diff --git a/CGParallel.h b/CGParallel.h
@@ -0,0 +1,22 @@
+/* 
+ * File:   CGParallel.h
+ * Author: zhakov
+ *
+ * Created on 13 Декабрь 2012 г., 21:37
+ */
+
+#ifndef CGPARALLEL_H
+#define	CGPARALLEL_H
+
+class CGParallel {
+public:
+    CGParallel();
+    CGParallel(const CGParallel& orig);
+    virtual ~CGParallel();
+    void resultCalculation(double** pMatrix, double* pVector, double* pResult, int Size);
+private:
+    double diff(double *vector1, double* vector2, int Size);
+};
+
+#endif	/* CGPARALLEL_H */
+
diff --git a/gaussParallel.cpp b/gaussParallel.cpp
@@ -32,9 +32,9 @@ int gaussParallel::resultCalculation(double** pMatrix, double* pVector, double*
 
 
     // Gaussian elimination
-    serialGaussianElimination(pMatrix, pVector);
+    gaussianElimination(pMatrix, pVector);
     // Back substitution
-    serialBackSubstitution(pMatrix, pVector, pResult);
+    backSubstitution(pMatrix, pVector, pResult);
 
     return 0;
 }
@@ -61,18 +61,25 @@ int gaussParallel::findPivotRow(double** pMatrix, int Iter) {
                 ThreadPivotRow.MaxValue = fabs(pMatrix[i][Iter]);
             }
         }
-        printf("Local thread(id = %i) pivot row: %i\n", omp_get_thread_num(), ThreadPivotRow.PivotRow);
-    }
+#pragma omp critical
+        {
+            if (ThreadPivotRow.MaxValue > MaxValue) {
+                MaxValue = ThreadPivotRow.MaxValue;
+                PivotRow = ThreadPivotRow.PivotRow;
+            }
+        } // pragma omp critical
+    }// pragma omp parallel
 
     return PivotRow;
 }
 
 
 // Column elimination
 
-int gaussParallel::serialColumnElimination(double** pMatrix, double* pVector, int Pivot, int Iter) {
+int gaussParallel::columnElimination(double** pMatrix, double* pVector, int Pivot, int Iter) {
     double PivotValue, PivotFactor;
     PivotValue = pMatrix[Pivot][Iter];
+#pragma omp parallel for private (PivotFactor)
     for (int i = 0; i < mSize; i++) {
         if (pSerialPivotIter[i] == -1) {
             PivotFactor = pMatrix[i][Iter] / PivotValue;
@@ -89,7 +96,7 @@ int gaussParallel::serialColumnElimination(double** pMatrix, double* pVector, in
 
 // Gaussian elimination
 
-int gaussParallel::serialGaussianElimination(double** pMatrix, double* pVector) {
+int gaussParallel::gaussianElimination(double** pMatrix, double* pVector) {
     int Iter;
     // The Number of the iteration of the gaussian
     // elimination
@@ -101,21 +108,22 @@ int gaussParallel::serialGaussianElimination(double** pMatrix, double* pVector)
         PivotRow = findPivotRow(pMatrix, Iter);
         pSerialPivotPos[Iter] = PivotRow;
         pSerialPivotIter[PivotRow] = Iter;
-        serialColumnElimination(pMatrix, pVector, PivotRow, Iter);
+        columnElimination(pMatrix, pVector, PivotRow, Iter);
     }
 
 
 
     return 0;
 }
 
-//Обратный ход метода Гаусса
-
-int gaussParallel::serialBackSubstitution(double** pMatrix, double* pVector, double* pResult) {
+/* Обратный ход метода Гаусса
+ */
+int gaussParallel::backSubstitution(double** pMatrix, double* pVector, double* pResult) {
     int RowIndex, Row;
     for (int i = mSize - 1; i >= 0; i--) {
         RowIndex = pSerialPivotPos[i];
         pResult[i] = pVector[RowIndex] / pMatrix[RowIndex][i];
+#pragma omp parallel for private (Row)
         for (int j = 0; j < i; j++) {
             Row = pSerialPivotPos[j];
             pVector[j] -= pMatrix[Row][i] * pResult[i];
diff --git a/gaussParallel.h b/gaussParallel.h
@@ -21,9 +21,9 @@ class gaussParallel {
     int* pSerialPivotIter; // The iterations, at which the rows were pivots
     
     int findPivotRow(double** pMatrix, int Iter);
-    int serialGaussianElimination(double** pMatrix, double* pVector);
-    int serialBackSubstitution(double** pMatrix, double* pVector, double* pResult);
-    int serialColumnElimination(double** pMatrix, double* pVector, int Pivot, int Iter);
+    int gaussianElimination(double** pMatrix, double* pVector);
+    int backSubstitution(double** pMatrix, double* pVector, double* pResult);
+    int columnElimination(double** pMatrix, double* pVector, int Pivot, int Iter);
 };
 
 #endif	/* GAUSSPARALLEL_H */
diff --git a/main.cpp b/main.cpp
@@ -126,7 +126,9 @@ int main(int argc, char** argv) {
 
 
 
-
+    //Проверяем результат
+    matrixHelpers::testSolvingResult(pMatrix, pVector, pResult, mSize);
+    
 
 
     //Потраченное время
diff --git a/matrixHelpers.cpp b/matrixHelpers.cpp
@@ -66,10 +66,10 @@ int matrixHelpers::testSolvingResult(double** pMatrix, double* pVector, double*
         }
     }
     if (equal == 1) {
-        printf("The result of the parallel Gauss algorithm is NOT correct."
-                "Check your code.");
+        printf("The result of the algorithm is NOT correct."
+                "Check your code.\n");
     } else {
-        printf("The result of the parallel Gauss algorithm is correct.");
+        printf("The result of the algorithm is correct.\n");
     }
     
     delete [] pRightPartVector;
diff --git a/nbproject/Makefile-Debug.mk b/nbproject/Makefile-Debug.mk
@@ -36,6 +36,7 @@ OBJECTDIR=${CND_BUILDDIR}/${CND_CONF}/${CND_PLATFORM}
 # Object Files
 OBJECTFILES= \
 	${OBJECTDIR}/main.o \
+	${OBJECTDIR}/CGParallel.o \
 	${OBJECTDIR}/matrixHelpers.o \
 	${OBJECTDIR}/dataGen.o \
 	${OBJECTDIR}/gaussSerial.o \
@@ -71,6 +72,11 @@ ${OBJECTDIR}/main.o: main.cpp
 	${RM} $@.d
 	$(COMPILE.cc) -g -MMD -MP -MF $@.d -o ${OBJECTDIR}/main.o main.cpp
 
+${OBJECTDIR}/CGParallel.o: CGParallel.cpp 
+	${MKDIR} -p ${OBJECTDIR}
+	${RM} $@.d
+	$(COMPILE.cc) -g -MMD -MP -MF $@.d -o ${OBJECTDIR}/CGParallel.o CGParallel.cpp
+
 ${OBJECTDIR}/matrixHelpers.o: matrixHelpers.cpp 
 	${MKDIR} -p ${OBJECTDIR}
 	${RM} $@.d
diff --git a/nbproject/Makefile-Release.mk b/nbproject/Makefile-Release.mk
@@ -36,6 +36,7 @@ OBJECTDIR=${CND_BUILDDIR}/${CND_CONF}/${CND_PLATFORM}
 # Object Files
 OBJECTFILES= \
 	${OBJECTDIR}/main.o \
+	${OBJECTDIR}/CGParallel.o \
 	${OBJECTDIR}/matrixHelpers.o \
 	${OBJECTDIR}/dataGen.o \
 	${OBJECTDIR}/gaussSerial.o \
@@ -71,6 +72,11 @@ ${OBJECTDIR}/main.o: main.cpp
 	${RM} $@.d
 	$(COMPILE.cc) -O2 -MMD -MP -MF $@.d -o ${OBJECTDIR}/main.o main.cpp
 
+${OBJECTDIR}/CGParallel.o: CGParallel.cpp 
+	${MKDIR} -p ${OBJECTDIR}
+	${RM} $@.d
+	$(COMPILE.cc) -O2 -MMD -MP -MF $@.d -o ${OBJECTDIR}/CGParallel.o CGParallel.cpp
+
 ${OBJECTDIR}/matrixHelpers.o: matrixHelpers.cpp 
 	${MKDIR} -p ${OBJECTDIR}
 	${RM} $@.d
diff --git a/nbproject/configurations.xml b/nbproject/configurations.xml
@@ -4,6 +4,8 @@
     <logicalFolder name="SourceFiles"
                    displayName="Исходные файлы"
                    projectFiles="true">
+      <itemPath>CGParallel.cpp</itemPath>
+      <itemPath>CGParallel.h</itemPath>
       <itemPath>dataGen.cpp</itemPath>
       <itemPath>dataGen.h</itemPath>
       <itemPath>gaussParallel.cpp</itemPath>
diff --git a/nbproject/private/configurations.xml b/nbproject/private/configurations.xml
@@ -24,18 +24,18 @@
       </nativedebugger>
       <runprofile version="9">
         <runcommandpicklist>
-          <runcommandpicklistitem>"${OUTPUT_PATH}"</runcommandpicklistitem>
           <runcommandpicklistitem>"${OUTPUT_PATH}" -m 64</runcommandpicklistitem>
           <runcommandpicklistitem>"${OUTPUT_PATH}" -m asdf</runcommandpicklistitem>
           <runcommandpicklistitem>"${OUTPUT_PATH}" -m 5</runcommandpicklistitem>
           <runcommandpicklistitem>"${OUTPUT_PATH}" -m 3000</runcommandpicklistitem>
           <runcommandpicklistitem>"${OUTPUT_PATH}" -m 500</runcommandpicklistitem>
           <runcommandpicklistitem>"${OUTPUT_PATH}" -m 500 -algorithm </runcommandpicklistitem>
           <runcommandpicklistitem>"${OUTPUT_PATH}" -m 500 -algorithm gauss-serial</runcommandpicklistitem>
-          <runcommandpicklistitem>"${OUTPUT_PATH}" -msize 500 -algorithm gauss-serial</runcommandpicklistitem>
           <runcommandpicklistitem>"${OUTPUT_PATH}" -msize 500 -algorithm gauss-parallel</runcommandpicklistitem>
+          <runcommandpicklistitem>"${OUTPUT_PATH}" -msize 500 -algorithm gauss-serial</runcommandpicklistitem>
+          <runcommandpicklistitem>"${OUTPUT_PATH}" -msize 500 -algorithm CG-parallel</runcommandpicklistitem>
         </runcommandpicklist>
-        <runcommand>"${OUTPUT_PATH}" -msize 500 -algorithm gauss-parallel</runcommand>
+        <runcommand>"${OUTPUT_PATH}" -msize 500 -algorithm CG-parallel</runcommand>
         <rundir></rundir>
         <buildfirst>true</buildfirst>
         <terminal-type>0</terminal-type>

Original file line number	Diff line number	Diff line change
`@@ -66,10 +66,10 @@ int matrixHelpers::testSolvingResult(double** pMatrix, double* pVector, double*`
`66`	`66`	`}`
`67`	`67`	`}`
`68`	`68`	`if (equal == 1) {`
`69`		`- printf("The result of the parallel Gauss algorithm is NOT correct."`
`70`		`- "Check your code.");`
	`69`	`+ printf("The result of the algorithm is NOT correct."`
	`70`	`+ "Check your code.\n");`
`71`	`71`	`} else {`
`72`		`- printf("The result of the parallel Gauss algorithm is correct.");`
	`72`	`+ printf("The result of the algorithm is correct.\n");`
`73`	`73`	`}`
`74`	`74`
`75`	`75`	`delete [] pRightPartVector;`