diff --git a/AMDiS/src/UmfPackSolver.hh b/AMDiS/src/UmfPackSolver.hh index c8007a7bd3420f5426c8ecc1fafe5472ccb171c1..9762b0fa10cd1ea688ed70bd5718ef32a95c3700 100644 --- a/AMDiS/src/UmfPackSolver.hh +++ b/AMDiS/src/UmfPackSolver.hh @@ -44,8 +44,22 @@ namespace AMDiS { // Number of systems. int nComponents = m->getSize(); - // Size of the new composed matrix. - int newMatrixSize = ((*m)[0][0])->getFESpace()->getAdmin()->getUsedSize() * nComponents; + + // Calculate size of the new composed matrix. + int newMatrixSize = 0; + std::vector<int> matricesSize(nComponents); + + for (int i = 0; i < nComponents; i++) { + for (int j = 0; j < nComponents; j++) { + if ((*m)[i][j]) { + matricesSize[i] = ((*m)[i][j])->getFESpace()->getAdmin()->getUsedSize(); + newMatrixSize += matricesSize[i]; + break; + } + } + } + + // int newMatrixSize = ((*m)[0][0])->getFESpace()->getAdmin()->getUsedSize() * nComponents; // The new matrix has to be stored in compressed col format, therefore // the cols are collected. @@ -54,34 +68,44 @@ namespace AMDiS { // Counter for the number of non-zero elements in the new matrix. int nElements = 0; - for (int stencilRow = 0; stencilRow < nComponents; stencilRow++) { + for (int stencilRow = 0, shiftRow = 0; stencilRow < nComponents; stencilRow++) { + + int shiftCol = 0; + for (int stencilCol = 0; stencilCol < nComponents; stencilCol++) { if (!(*m)[stencilRow][stencilCol]) { + shiftCol += matricesSize[stencilCol]; continue; } - + DOFMatrix::Iterator matrixRow((*m)[stencilRow][stencilCol], USED_DOFS); int rowIndex = 0; for (matrixRow.reset(); !matrixRow.end(); matrixRow++, rowIndex++) { - for (int i = 0; i < static_cast<int>((*matrixRow).size()); i++) { + for (int i = 0; i < static_cast<int>((*matrixRow).size()); i++) { if ((*matrixRow)[i].col >= 0) { MatEntry me; me.entry = (*matrixRow)[i].entry; // The col field is used to store the row number of the new element. - me.col = (rowIndex * nComponents) + stencilRow; + me.col = rowIndex + shiftRow; // And save the new element in the corresponding column. - cols[((*matrixRow)[i].col * nComponents) + stencilCol].push_back(me); + cols[(*matrixRow)[i].col + shiftCol].push_back(me); nElements++; } - } - } + } + } + + + shiftCol += matricesSize[stencilCol]; } + + shiftRow += matricesSize[stencilRow]; } + // Data fields for UMFPack. int *Ap = (int*)malloc(sizeof(int) * (newMatrixSize + 1)); int *Ai = (int*)malloc(sizeof(int) * nElements); @@ -90,15 +114,15 @@ namespace AMDiS { double *xvec = (double*)malloc(sizeof(double) * newMatrixSize); // Resort the right hand side of the linear system. - for (int i = 0; i < b->getSize(); i++) { + for (int i = 0, counter = 0; i < b->getSize(); i++) { DOFVector<double>::Iterator it(b->getDOFVector(i), USED_DOFS); - int counter = 0; - for (it.reset(); !it.end(); ++it, counter++) { - bvec[counter * nComponents + i] = *it; + for (it.reset(); !it.end(); ++it) { + bvec[counter++] = *it; } } + // Create fields Ap, Ai and Ax. int elCounter = 0; Ap[0] = 0; @@ -164,12 +188,11 @@ namespace AMDiS { // Copy and resort solution. - for (int i = 0; i < x->getSize(); i++) { + for (int i = 0, counter = 0; i < x->getSize(); i++) { DOFVector<double>::Iterator it(x->getDOFVector(i), USED_DOFS); - - int counter = 0; - for (it.reset(); !it.end(); it++, counter++) { - *it = xvec[counter * nComponents + i]; + + for (it.reset(); !it.end(); it++) { + *it = xvec[counter++]; } }