Added SolveTask, LUfactorization(), device-memory-specific code

Include/dpsim/MNASolverGpu.h

@@ -25,22 +25,44 @@ namespace DPsim {
 		cusolverDnHandle_t mCusolverHandle;
 		///Stream
 		cudaStream_t mStream;
-		/// Device copy of System-Matrix
-		double *mGpuSystemMatrix;
-		/// Device copy of Right Vector
-		double *mGpuRightVector;
-		/// Device copy if Left Vector
-		double *mGpuLeftVector;
+
+		/// Variables for solving one Equation-system
+		struct GpuData {
+			/// Device copy of System-Matrix
+			double *matrix;
+			/// Device copy of Vector
+			double *rightVector;
+
+			/// Device-Workspace for getrf
+			double *workSpace;
+			/// Pivoting-Sequence
+			int *pivSeq;
+			/// Errorinfo
+			int *errInfo;
+		} mDeviceCopy;
 
 		/// Initialize cuSolver-library
         void initialize();
         /// Allocate Space for Vectors & Matrices on GPU
-        void createEmptyVectors();
-	    void createEmptySystemMatrix();
+        void allocateDeviceMemory();
+		/// Copy Systemmatrix to Device
+		void copySystemMatrixToDevice();
+		/// LU factorization
+		void LUfactorization();
+		
 
 	public:
 		MnaSolverGpu();
 
 		virtual ~MnaSolverGpu();
+
+
+		class SolveTask : public MnaSolver::SolveTask {
+			SolveTask(MnaSolver<VarType>& solver, Bool steadyStateInit) :
+			MnaSolver::SolveTask(MnaSolver<VarType>& solver, Bool steadyStateInit) {
+			}
+
+			void execute(Real time, Int timeStepCount);
+		};
     };
 }
\ No newline at end of file

Source/MNASolverGpu.cpp

@@ -7,9 +7,9 @@ namespace DPsim {
 
 template <typename VarType>
 MnaSolverGpu<VarType>::MnaSolverGpu() :
-mCusolverHandle(nullptr), mStream(nullptr),
-mGpuSystemMatrix(nullptr), mGpuLeftVector(nullptr), mGpuRightVector(nullptr) {
+mCusolverHandle(nullptr), mStream(nullptr) {
 
+    mDeviceCopy = {nullptr, nullptr, nullptr, nullptr, nullptr, nullptr};
     cusolverStatus_t status = CUSOLVER_STATUS_SUCCESS;
     cudaError_t cudaErrorCode;
 
@@ -27,13 +27,17 @@ MnaSolverGpu<VarType>::~MnaSolverGpu() {
     if(mStream)
         cudaStreamDestroy(stream);
 
-    //Matrix & Vectors
-    if(mGpuSystemMatrix)
-        cudaFree(mGpuSystemMatrix);
-    if(mGpuLeftVector)
-        cudaFree(mGpuLeftVector);
-    if(mGpuRightVector)
-        cudaFree(mGpuRightVector);
+    //Memory allocated on device
+    if(mDeviceCopy.matrix)
+        cudaFree(mDeviceCopy.matrix);
+    if(mDeviceCopy.rightVector)
+        cudaFree(mDeviceCopy.rightVector);
+    if(mDeviceCopy.workSpace)
+        cudaFree(mDeviceCopy.workSpace);
+    if(mDeviceCopy.pivSeq)
+        cudaFree(mDeviceCopy.pivSeq);
+    if(mDeviceCopy.errInfo)
+        cudaFree(mDevice.errInfo);
 
     cudaDeviceReset();
 }
@@ -42,27 +46,113 @@ template <typename VarType>
 void MnaSolverGpu<VarType>::initialize() {
     MnaSolver::initialize();
 
-    createEmptyVectors();
-    createEmptySystemMatrix();
-
+    allocateDeviceMemory();
     //Copy Systemmatrix to device
+    copySystemMatrixToDevice();
+    //LU factorization
 }
 
 /// Allocate Space for Vectors & Matrices
 template <typename VarType>
-void MnaSolverGpu<VarType>::createEmptyVectors() {
-    //TODO Error Checking
+void MnaSolverGpu<VarType>::allocateDeviceMemory() {
+    //TODO Error checking
+    //Vectors
     auto size = sizeof(Real) * mNumSimNodes;
     cudaError_t stat;
-    stat = cudaMalloc(static_cast<void**>(&mGpuLeftVector), size);
-    stat = cudaMalloc(static_cast<void**>(&mGpuRightVector), size);
+    stat = cudaMalloc(static_cast<void**>(&mDeviceCopy.leftVector), size);
+    stat = cudaMalloc(static_cast<void**>(&mDeviceCopy.rightVector), size);
+    //Matrix
+    size = sizeof(Real) * mNumSimNodes * mNumSimNodes;
+    stat = cudaMalloc(static_cast<void**>(&mDeviceCopy.matrix), size);
+    //Pivoting-Sequence
+    stat = cudaMalloc(static_cast<void**>(&mDeviceCopy.matrix), sizeof(int) * mNumSimNodes);
+    //Errorcode
+    stat = cudaMalloc(static_cast<void**>(&mDeviceCopy.errInfo), sizeof(int));
+
+    //Workspace
+    int workSpaceSize = 0;
+    cusolverStatus_t status = cusolverDnDgetrf_bufferSize(
+        mCusolverHandle,
+        mNumSimNodes,
+        mNumSimNodes,
+        mDeviceCopy.matrix,
+        mNumSimNodes,
+        &workSpaceSize);
+    stat = cudaMalloc(&mDeviceCopy.workSpace, workSpaceSize);
 }
 
 template <typename VarType>
-void MnaSolverGpu<VarType>::createEmptySystemMatrix() {
+void MnaSolverGpu<VarType>::copySystemMatrixToDevice() {
     //TODO Error Checking
-    auto size = sizeof(Real) * mNumSimNodes * mNumSimNodes;
-    cudaError_t stat = cudaMalloc(static_cast<void**>(&mGpuLeftVector), size);
+    Real *mat = new Real[mNumSimNodes * mNumSimNodes];
+    for(int i = 0; i < mNumSimNodes; i++) {
+        for(int j = 0; j < mNumSimNodes; j++) {
+            mat[i * mNumSimNodes + j] = systemMatrix()(i, j);
+        }
+    }
+    cudaError_t cudaStat1 = cudaMemcpy(mDeviceCopy.matrix, mat, sizeof(Real) * mNumSimNodes * mNumSimNodes, cudaMemcpyHostToDevice);
+}
+
+template <typename VarType>
+void MnaSolverGpu<VarType>::LUfactorization() {
+    //TODO Error checking
+    auto status = cusolverDnDgetrf(
+        mCusolverHandle,
+        mNumSimNodes,
+        mNumSimNodes,
+        mDeviceCopy.matrix,
+        mNumSimNodes,
+        mDeviceCopy.workSpace,
+        mDeviceCopy.pivSeq,
+        mDeviceCopy.errInfo);
+    status = cudaDeviceSynchronize();
+}
+
+template <typename VarType>
+void MnaSolverGpu<VarType>::SolveTask::execute(Real time, Int timeStepCount) {
+    // Reset source vector
+	mSolver.mRightSideVector.setZero();
+
+    // Add together the right side vector (computed by the components'
+	// pre-step tasks)
+	for (auto stamp : mSolver.mRightVectorStamps)
+		mSolver.mRightSideVector += *stamp;
+    //Copy to device
+    double *buffer = new double[mNumSimNodes];
+    for(int i = 0; i < mNumSimNodes; i++) {
+        buffer[i] = mSolver.mRightSideVector(1, i); //TODO check
+    }
+    cudaError_t ec = cudaMemcpy(mDeviceCopy.rightVector, buffer, mNumSimNodes, cudaMemcpyHostToDevice);
+
+    // Solve
+	if (mSolver.mSwitchedMatrices.size() > 0) {
+        auto status = cusolverDnDgetrs(
+            mCusolverHandle,
+            CUBLAS_OP_N,
+            mNumSimNodes,
+            1, /* nrhs */
+            mDeviceCopy.matrix,
+            mNumSimNodes,
+            mDeviceCopy.pivSeq,
+            mDeviceCopy.rightVector,
+            mNumSimNodes,
+            mDeviceCopy.errInfo);
+    }
+
+    //Copy Leftvector back
+    ec = cudaMemcpy(buffer, mDeviceCopy.rightVector, mNumSimNodes, cudaMemcpyDeviceToHost);
+    for(int i = 0; i < mNumSimNodes; i++) {
+        mSolver.mLeftSideVector(1, i) = buffer[i]; // TODO check
+    }
+
+	// TODO split into separate task? (dependent on x, updating all v attributes)
+	for (UInt nodeIdx = 0; nodeIdx < mSolver.mNumNetNodes; nodeIdx++)
+		mSolver.mNodes[nodeIdx]->mnaUpdateVoltage(mSolver.mLeftSideVector);
+
+	if (!mSteadyStateInit)
+		mSolver.updateSwitchStatus();
+
+	// Components' states will be updated by the post-step tasks
 }
 
 }