System of equation solver

system/main.cc

+// ------------------------------------------------------------
+// main.cc
+//
+// ------------------------------------------------------------
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <complex>
+#include <iostream>
+#include <string>
+#include <vector>
+#include "../program2D/source/ugblock2D.h"
+#include "source/vectorVar.h"
+#include "solver/iterativeSolversSystem.h"
+
+
+using std::complex;
+using namespace ::_COLSAMM_;
+
+#define L2NORM(vector) sqrt(product(vector,vector).real())
+#define RESTART 0
+
+
+double isZero(double x)
+{
+	return fabs(x) < 1e-10?1.0:0.0;
+}
+
+
+
+int main(int argc, char** argv) {
+	std::ofstream DATEI;
+
+	int n = 10;
+
+	//SectionedRectangle geo(v_x, v_y);  
+        //Rectangle geo(1.0, 1.0);
+	Disc geo(2.0,1.0);
+
+	Boundary_Marker2D boundary(&geo);
+	Unstructured2DGrid_Marker2D inner_points(&geo);
+	inner_points.complement(&boundary);
+
+        // boudary integral test 
+        Blockgrid2D  blockBoTest(&geo,n);
+        Variable2D<double> one(blockBoTest);
+        Variable2D<double>  ff(blockBoTest);
+
+        Variable2D<double>  u(blockBoTest);
+        Variable2D<double>  interpol(blockBoTest);
+        Variable2D<double>  g(blockBoTest);
+        Variable2D<double>  p(blockBoTest);	
+	
+        Local_stiffness_matrix2D<double> boIntegral(blockBoTest);  
+
+	DATEI.open("ug.vtk",std::ios :: out);
+	geo.PrintVtk(&DATEI);
+	DATEI.close();
+
+	X_coordinate2d X(blockBoTest);
+	Y_coordinate2d Y(blockBoTest);
+
+/*
+        Local_stiffness_matrix2D<double>  locStiff(blockBoTest);
+	//	Test.Calculate(v_()*w_());
+        locStiff.Calculate_HelmBoundary();        
+	one = 1.0;
+	ff = locStiff(one);	
+	cout << "inegral Rand: " << product(ff,one) << endl;
+*/
+
+        u = X;
+        g = Y -0.2;    
+
+        VariableVector<double, Variable2D<double> > vectorU(u,u);
+        VariableVector<double, Variable2D<double> > vectorB(g,g);
+        VariableVector<double, Variable2D<double> > vectorI(interpol,p);  
+        vectorI = vectorU + vectorB;
+  
+        cout << " Maximum I: " << Maximum(vectorI) << endl;
+        cout << " Minimum I: " << Minimum(vectorI) << endl;
+
+        cout << " Maximum U: " << Maximum(vectorU) << endl;
+        cout << " Minimum U: " << Minimum(vectorU) << endl;
+
+        cout << " Maximum B: " << Maximum(vectorB) << endl;
+        cout << " Minimum B: " << Minimum(vectorB) << endl;
+	
+	
+	DATEI.open("interpol.vtk");
+	interpol.Print_VTK(DATEI);
+	DATEI.close();	
+
+	DATEI.open("one.vtk");
+	one.Print_VTK(DATEI);
+	DATEI.close();
+
+	DATEI.open("ff.vtk");
+	ff.Print_VTK(DATEI);
+	DATEI.close();
+}
+

system/para.dat

+2 5 2

system/program.pro

+######################################################################
+# Automatically generated by qmake (2.01a) Di. Sep 16 15:38:39 2008
+######################################################################
+
+# getting pathes:
+include(../config.pri)
+
+# GPU support
+#CUDA_BASE = /usr/local
+#GPUSOLVE_PATH = ../cuda/GPUsolve-0.1.0/
+
+#QMAKE_CXXFLAGS  += -I/usr/local/include -I/usr/local/cuda/include/
+#DEFINES += _GPU 
+#LIBS += -Lsource/cuda/lib -L/usr/local/cuda/lib/ -lcuda -lcudart -lgpusolve
+
+
+# The code is too messy: lots of unused-paramter warnings,
+# so turned off here:
+unix:QMAKE_CXXFLAGS_WARN_ON += -Wno-unused-parameter
+#win32:QMAKE_CXXFLAGS_WARN_ON += /w44100
+win32:QMAKE_CXXFLAGS_WARN_ON = /W3 /w44100 /w44503
+#unix:QMAKE_CXXFLAGS+="-fexec-charset=ISO-8859-1"
+
+#CONFIG -= openmp
+#CONFIG += openmp
+CONFIG += debug 
+
+# placing *.o files in bin_o
+unix: OBJECTS_DIR= bin_o
+unix: UI_DIR =  UI_FILES
+unix: MOC_DIR =  MOC_FILES
+
+
+# OpenMP support
+
+openmp {
+    # gcc
+    unix:QMAKE_CXXFLAGS  += -fopenmp
+    # Windows with msvc cl
+    win32:QMAKE_CXXFLAGS  += /openmp
+    LIBS += $$OPENMP_LIBPATH $$OPENMP_LIBS
+} else  {
+    # prevent warnings about unknown omp pragmas
+    # WARNING: This also disables warnings about other pragmas.
+#   unix:QMAKE_CXXFLAGS_WARN_ON+=-Wno-unknown-pragmas
+}
+
+				
+
+include(system.pri)
+include(../program/ugblocks3D.pri)
+include(../program2D/ugblocks2D.pri)
+include(../colsamm.pri)
+
+
+SOURCES += main.cc

system/solver/gmres/gmresSystem.h

+/**********************************************************************************
+ * Copyright 2020 Christoph Pflaum
+ * 		Department Informatik Lehrstuhl 10 - Systemsimulation
+ *		Friedrich-Alexander Universität Erlangen-Nürnberg
+ * 
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ **********************************************************************************/
+
+//#define _USE_MATH_DEFINES
+#include <cmath>
+#include <complex>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+using std::vector;
+
+
+#ifndef GMRESSYSTEM_H
+#define GMRESSYSTEM_H
+
+/**
+ GMRES Löser für \f$ A x = b \f$
+ **/
+template<class DTyp, class Vari, class Mat>
+  class GmResSystem : public IterativeSolverSystemAbstract <DTyp,Vari,Mat> {
+
+  public:
+    /** 
+     * @param blockgrid Gitter auf dem die Gleichung gelösst werden soll
+     **/
+    GmResSystem(VariableVector<DTyp,Vari>* res);
+    //~IterativeSolverAbstract() { };
+    ~GmResSystem();
+
+    /** 
+     * set maximum size of Kryloc-space
+     * @param numberRestart_ important GMRES parameter, which is default 100, changing this parameter changes memory allocation
+     **/
+    void setRestartNumber(int numberRestart_);
+    
+    /**
+     * solves the linear equation system \f$ A x = b \f$
+     * class Mat must contain member function   Mat::apply(x,b)
+     **/  
+    void solve(Mat& A, VariableVector<DTyp,Vari>& b, VariableVector<DTyp,Vari>& x);
+
+    bool isResiduumCalculated() { return true; }    
+  private:
+
+    std::vector< VariableVector<DTyp,Vari>* > Q;
+
+    int Restart;
+
+    DTyp** H;     // size [Restart + 1][Restart]
+    DTyp* C;      // size [Restart]
+    DTyp* S;      // size [Restart]
+    DTyp* y;      // size [Restart]    
+    DTyp* gamma;  // size [Restart+1]
+
+    static void GivensRotation(const complex<double>& x1, complex<double>& x2, complex<double>& c, complex<double>& s);    
+    static void GivensRotation(const double& x1, double& x2, double& c, double& s);        
+};
+
+
+ 
+  
+#endif

system/solver/gmres/gmresSystem_cc.h

+
+/**********************************************************************************
+ * Copyright 2020 Christoph Pflaum
+ * 		Department Informatik Lehrstuhl 10 - Systemsimulation
+ *		Friedrich-Alexander Universität Erlangen-Nürnberg
+ * 
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ **********************************************************************************/
+
+
+
+template<class DTyp, class Vari, class Mat>
+GmResSystem<DTyp,Vari,Mat>::GmResSystem(VariableVector<DTyp,Vari>* res) : IterativeSolverSystemAbstract<DTyp,Vari,Mat>(res) {
+   Restart = 100;
+   
+   Q.resize(Restart+1);
+   
+   H = new DTyp*[Restart+1];
+   C = new DTyp[Restart];
+   S = new DTyp[Restart];
+   y = new DTyp[Restart];
+   gamma = new DTyp[Restart+1];
+   
+   for(int i=0;i<Restart+1;++i) {
+       H[i] = new DTyp[Restart];
+       if(this->size==1) {
+	  Vari* v = new Vari(*this->rXcomponent);
+	  Q.at(i) = new VariableVector<DTyp,Vari>(*v);
+       }
+       if(this->size==2) {
+	  Vari* vx = new Vari(*this->rXcomponent);
+	  Vari* vy = new Vari(*this->rXcomponent);	 
+	  Q.at(i) = new VariableVector<DTyp,Vari>(*vx,*vy);
+       }
+       if(this->size==3) {
+	  Vari* vx = new Vari(*this->rXcomponent);
+	  Vari* vy = new Vari(*this->rXcomponent); 
+	  Vari* vz = new Vari(*this->rXcomponent);	  
+	  Q.at(i) = new VariableVector<DTyp,Vari>(*vx,*vy,*vz);
+       }
+   }
+}
+
+template<class DTyp, class Vari, class Mat>
+GmResSystem<DTyp,Vari,Mat>::~GmResSystem() {
+   for(int i=0;i<Restart+1;++i) {
+       delete[] H[i];
+       Q.at(i)->deleteVariables();
+       delete Q.at(i);
+   }
+   delete H;
+   delete C;
+   delete S;
+   delete y;
+   delete gamma;
+} 
+
+
+template<class DTyp, class Vari, class Mat>
+void GmResSystem<DTyp,Vari,Mat>::setRestartNumber(int numberRestart_) {
+  int oldRestart = Restart;
+  Restart = numberRestart_;
+  if(oldRestart==Restart) return;
+
+  for(int i=0;i<oldRestart+1;++i) delete H[i];
+  delete H;
+  delete C;
+  delete S;
+  delete y;
+  delete gamma;
+    
+   H = new DTyp*[Restart+1];
+   C = new DTyp[Restart];
+   S = new DTyp[Restart];
+   y = new DTyp[Restart];
+   gamma = new DTyp[Restart+1];
+   
+   for(int i=0;i<Restart+1;++i) {
+       H[i] = new DTyp[Restart];
+   }
+  
+  
+  if(oldRestart<Restart) {
+     Q.resize(Restart+1);
+     for(int i=oldRestart+1;i<Restart+1;++i) {
+         //Q.at(i) = new Variable2D<DTyp>(*(this->blockgrid));
+         if(this->size==1) {
+	    Vari* v = new Vari(*this->rXcomponent);
+	    Q.at(i) = new VariableVector<DTyp,Vari>(*v);
+         }
+         if(this->size==2) {
+	    Vari* vx = new Vari(*this->rXcomponent);
+	    Vari* vy = new Vari(*this->rXcomponent);	 
+	    Q.at(i) = new VariableVector<DTyp,Vari>(*vx,*vy);
+         }
+         if(this->size==3) {
+	    Vari* vx = new Vari(*this->rXcomponent);
+	    Vari* vy = new Vari(*this->rXcomponent); 
+	    Vari* vz = new Vari(*this->rXcomponent);	  
+	    Q.at(i) = new VariableVector<DTyp,Vari>(*vx,*vy,*vz);
+         }	 
+     }
+  }
+  else {
+     for(int i=Restart+1;i<oldRestart+1;++i) {
+         Q.at(i)->deleteVariables();
+         delete Q.at(i);
+     }
+     Q.resize(Restart+1);
+  }
+}
+
+
+
+template<class DTyp, class Vari, class Mat>
+inline void GmResSystem<DTyp, Vari, Mat>::GivensRotation(const complex<double>& x1, complex<double>& x2, 
+						    complex<double>& c, complex<double>& s)  {
+    complex<double> rho = x1 / abs(x1) * sqrt(x1.real() * x1.real() + x1.imag() * x1.imag() + x2.real() * x2.real() + x2.imag() * x2.imag());
+    c = x1 / rho;
+    s = std::conj(x2 / rho);
+  }
+
+template<class DTyp, class Vari, class Mat>
+inline void GmResSystem<DTyp, Vari, Mat>::GivensRotation(const double& x1, double& x2, double& c, double& s)  {
+    double rho = x1 / fabs(x1) * sqrt(x1 * x1 + x2 * x2);
+    c = x1 / rho;
+    s = x2 / rho;
+  }
+
+  
+
+
+
+template<class DTyp, class Vari, class Mat>
+void GmResSystem<DTyp,Vari,Mat>::solve(Mat& A, VariableVector<DTyp,Vari>& b, VariableVector<DTyp,Vari>& x) {
+    this->numberIterationsPerformed = 0;
+      
+//test
+//cout << " max_iter: " << this->max_iter << endl;    
+    
+    double r_abs;
+    do {
+      (*this->r) = 0.0;
+      A.apply(x,*this->r);
+      (*this->r) = b - (*this->r);      //  r = (b - A(x));
+     
+      r_abs = L2Norm(*this->r);
+      if(r_abs < 1e-12) break;
+
+      gamma[0] = r_abs;
+      *Q[0] = (*this->r) / gamma[0];
+
+
+      
+      for(int j = 1; j <= Restart; ++j) {
+        this->numberIterationsPerformed++;
+        *Q[j] = 0.0;
+	A.apply(*Q[j-1] , *Q[j]);  // *Q[j] = (A(*Q[j-1]));
+
+        for (int i = 0; i < j; ++i) {
+          H[i][j - 1] = product(*Q[j], *Q[i]);
+          *Q[j] = *Q[j] - H[i][j - 1] * *Q[i];
+        }
+
+        //H[j][j - 1] = L2NORM(*Q[j]);  // old
+        H[j][j - 1] = L2Norm(*Q[j]);    
+        *Q[j] = *Q[j] / H[j][j - 1];
+
+        for (int i = 1; i < j; ++i) {
+          const DTyp h_tmp = H[i - 1][j - 1];
+          H[i - 1][j - 1] = (C[i - 1]) * h_tmp + S[i - 1] * H[i][j - 1];
+          H[i][j - 1] = -std::conj(S[i - 1]) * h_tmp + C[i - 1] * H[i][j - 1];
+        }
+      
+        GivensRotation(H[j - 1][j - 1], H[j][j - 1], C[j - 1], S[j - 1]);
+
+        H[j - 1][j - 1] = (C[j - 1]) * H[j - 1][j - 1] + S[j - 1] * H[j][j - 1];
+
+        gamma[j] = -std::conj(S[j - 1]) * gamma[j - 1];
+        gamma[j - 1] = (C[j - 1]) * gamma[j - 1];
+        r_abs = abs(gamma[j]);
+
+        if (j == Restart || r_abs < this->eps || this->numberIterationsPerformed >= this->max_iter) {
+          for (int i = j - 1; i >= 0; i--) {
+            DTyp tmp_sum = gamma[i];
+            for (int l = i + 1; l < j; ++l) {
+              tmp_sum -= H[i][l] * y[l];
+
+            }
+            y[i] = tmp_sum / H[i][i];
+
+            x = x + y[i] * *Q[i];
+          }
+          if (r_abs < this->eps || this->numberIterationsPerformed >= this->max_iter) {
+	    this->residuum = r_abs;
+	    return;
+	  }
+        }
+      }
+      
+    } while (r_abs > this->eps && this->numberIterationsPerformed < this->max_iter);
+    this->residuum = r_abs;
+    return;
+}
+
+
+  
+  
+  
+  
\ No newline at end of file

system/solver/iterativeSolversSystem.h

+/**********************************************************************************
+ * Copyright 2020 Christoph Pflaum
+ * 		Department Informatik Lehrstuhl 10 - Systemsimulation
+ *		Friedrich-Alexander Universität Erlangen-Nürnberg
+ * 
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ **********************************************************************************/
+
+
+#ifndef ITERATIVESOLVERSYSTEM_H
+#define ITERATIVESOLVERSYSTEM_H
+
+
+/**
+ Abstrake Klasse zur Beschreibung eines iterativen Lösers zur Lösung von \f$ A x = b \f$
+ **/
+template<class DTyp, class Vari, class Mat>
+class IterativeSolverSystemAbstract  {
+   public:
+        /** 
+	 * @param res
+	 **/
+        IterativeSolverSystemAbstract(VariableVector<DTyp,Vari>* res)  {  
+	       eps = 0.01;    stopByEps = true;    max_iter = 50;
+	       numberIterationsPerformed = 0;
+	       r    = res;
+	       rXcomponent = r->getComponent(0);
+	       size = r->getSize();
+	}
+	virtual ~IterativeSolverSystemAbstract() {
+	  //test NOW
+//cout << " ~IterativeSolverAbstract() " << endl;
+	  
+	}
+	 
+	/**
+	 * @param stopByEps falls true, endet die iteration durch Abbruchkriterium mit eps
+	 * @param eps       Abbruchkriterium
+	 * @param max_iter  maximale Anzahl der Iterationen
+	 * */
+	void setBasicParameter(double eps_, bool stopByEps_ = true, int max_iter_ = 500) { eps = eps_; stopByEps = stopByEps_; max_iter = max_iter_; }
+		
+	/**
+	 * solves the linear equation system \f$ A x = b \f$
+	 * class Mat must contain member function   Mat::apply(x,b)
+	 **/
+	virtual void solve(Mat& A, VariableVector<DTyp,Vari>& b, VariableVector<DTyp,Vari>& x) = 0;
+	
+	
+	/**
+	 * @return number of iterations which were performed
+	 **/
+	int numberOfPerformedIterations() { return numberIterationsPerformed; }
+	
+	/**
+	 * @return residuum after iteration. This is not calculated by every iteration method
+	 **/
+	double getResiduum() { return residuum; }
+	
+         /**
+	 * @return true when iteration method calculates residuum
+	 **/
+	virtual bool isResiduumCalculated() = 0;
+	
+protected:
+     double eps;
+     bool stopByEps;
+     int max_iter;
+     
+     int numberIterationsPerformed;
+     double residuum;
+     
+     VariableVector<DTyp,Vari>* r;
+     Vari* rXcomponent;     
+     int size;
+};
+
+#include "gmres/gmresSystem.h"
+#include "gmres/gmresSystem_cc.h"
+
+#endif //ITERATIVESOLVERSYSTEM_H

system/source/extempSystem.h

+/**********************************************************************************
+ * Copyright 2020 Christoph Pflaum 
+ * 		Department Informatik Lehrstuhl 10 - Systemsimulation
+ *		Friedrich-Alexander Universität Erlangen-Nürnberg
+ * 
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ **********************************************************************************/
+
+
+#ifndef EXTEMPSYSTEM_H
+#define EXTEMPSYSTEM_H
+
+
+
+template <class A> struct ExprSystem {
+	inline operator const A&() const {
+		return *static_cast<const A*> ( this );
+	}
+};
+
+//----------------------------------------------------------------------------
+template <class A, class B>
+class AddSystem : public ExprSystem<AddSystem<A, B> > {
+		const A& a_;
+		const B& b_;
+		
+	public:
+		inline AddSystem ( const A& a, const B& b ) : a_ ( a ), b_ ( b ) {}
+		
+						
+		typedef _TypeOf2_ ( A, B ) Result;
+		
+		template <int numComponent>
+		inline Result Give_data(int i) const {
+			return a_.template Give_data<numComponent> (i) +
+			       b_.template Give_data<numComponent> (i);
+		}		
+};
+
+//----------------------------------------------------------------------------
+template <class A, class B>
+inline AddSystem<A, B> operator+ ( const ExprSystem<A>& a, const ExprSystem<B>& b ) {
+	return AddSystem<A, B> ( a, b );
+}
+
+
+