[Getdp] PETSC --with-64-bit-indices=1
Mark van Doesburg
mark.van.doesburg at technolution.nl
Fri Jun 5 10:05:05 CEST 2009
Hi,
I've created a patch to be able to use PETSC-3.0.0p5 compiled with the
option --with-64-bit-indices=1. I don't actually have >2^31 degrees of
freedom, but it turns out that many algorithms fail on large problems
if PETSC is not compiled with 64-bit indices.
I'll try to port the Parallelization patch by Davide Rondini to the
latest snapshot.
Mark.
-------------- next part --------------
--- orig/getdp-2.0.0-cvs-20090603/Legacy/LinAlg_PETSC.cpp 2009-03-18 00:02:33.000000000 +0100
+++ getdp-2.0.0-cvs-20090603/Legacy/LinAlg_PETSC.cpp 2009-06-04 17:40:11.000000000 +0200
@@ -13,6 +13,7 @@
#include "LinAlg.h"
#include "MallocUtils.h"
#include "Message.h"
+#include <cstring>
// Johan, we curse you for a thousand generations!
#include "ProData.h"
@@ -68,9 +69,15 @@
// get additional petsc options from specified file (useful e.g. on
// Windows where we don't know where to search for ~/.petsrc)
+#if PETSC_VERSION_MAJOR==3
+ for(int i = 0; i < *argc - 1; i++)
+ if (!strcmp((*argv)[i], "-solver"))
+ PetscOptionsInsertFile(PETSC_COMM_WORLD, (*argv)[i+1], PETSC_FALSE);
+#else
for(int i = 0; i < *argc - 1; i++)
if (!strcmp((*argv)[i], "-solver"))
PetscOptionsInsertFile((*argv)[i+1]);
+#endif
}
void LinAlg_Finalize()
@@ -112,7 +119,7 @@
// database (if any)
ierr = MatSetFromOptions(M->M); MYCHECK(ierr);
// preallocation option must be set after other options
- int prealloc = 200;
+ PetscInt prealloc = 200;
PetscTruth set;
PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set);
ierr = MatSeqAIJSetPreallocation(M->M, prealloc, PETSC_NULL); MYCHECK(ierr);
@@ -182,9 +189,9 @@
void LinAlg_ScanVector(FILE *file, gVector *V)
{
- int n;
+ PetscInt n;
ierr = VecGetSize(V->V, &n); MYCHECK(ierr);
- for(int i = 0; i < n; i++){
+ for(PetscInt i = 0; i < n; i++){
double a, b;
PetscScalar tmp;
#if PETSC_USE_COMPLEX
@@ -210,11 +217,11 @@
void LinAlg_ReadVector(FILE *file, gVector *V)
{
- int n;
+ PetscInt n;
ierr = VecGetSize(V->V, &n); MYCHECK(ierr);
PetscScalar *tmp = (PetscScalar*)Malloc(n*sizeof(PetscScalar));
fread(tmp, sizeof(PetscScalar), n, file);
- for(int i = 0; i < n; i++){
+ for(PetscInt i = 0; i < n; i++){
ierr = VecSetValues(V->V, 1, &i, &tmp[i], INSERT_VALUES); MYCHECK(ierr);
}
Free(tmp);
@@ -236,7 +243,7 @@
void LinAlg_PrintVector(FILE *file, gVector *V)
{
- int n;
+ PetscInt n;
ierr = VecGetLocalSize(V->V, &n); MYCHECK(ierr);
PetscScalar *tmp;
ierr = VecGetArray(V->V, &tmp); MYCHECK(ierr);
@@ -267,7 +274,7 @@
void LinAlg_WriteVector(FILE *file, gVector *V)
{
- int n;
+ PetscInt n;
ierr = VecGetLocalSize(V->V, &n); MYCHECK(ierr);
PetscScalar *tmp;
ierr = VecGetArray(V->V, &tmp); MYCHECK(ierr);
@@ -283,22 +290,41 @@
void LinAlg_GetVectorSize(gVector *V, int *i)
{
- ierr = VecGetSize(V->V, i); MYCHECK(ierr);
+ PetscInt t;
+ ierr = VecGetSize(V->V, &t); MYCHECK(ierr);
+ if(t>INT_MAX)
+ Msg::Error("Problem too big\n");
+ *i=t;
}
void LinAlg_GetLocalVectorRange(gVector *V, int *low, int *high)
{
- ierr = VecGetOwnershipRange(V->V, low, high); MYCHECK(ierr);
+ PetscInt tlow, thigh;
+ ierr = VecGetOwnershipRange(V->V, &tlow, &thigh); MYCHECK(ierr);
+ if(tlow>INT_MAX || thigh>INT_MAX)
+ Msg::Error("Problem too big");
+ *low=tlow;
+ *high=thigh;
}
void LinAlg_GetMatrixSize(gMatrix *M, int *i, int *j)
{
- ierr = MatGetSize(M->M, i, j); MYCHECK(ierr);
+ PetscInt ti, tj;
+ ierr = MatGetSize(M->M, &ti, &tj); MYCHECK(ierr);
+ if(ti>INT_MAX || tj>INT_MAX)
+ Msg::Error("Problem too big");
+ *i=ti;
+ *j=tj;
}
void LinAlg_GetLocalMatrixRange(gMatrix *M, int *low, int *high)
{
- ierr = MatGetOwnershipRange(M->M, low, high); MYCHECK(ierr);
+ PetscInt tlow, thigh;
+ ierr = MatGetOwnershipRange(M->M, &tlow, &thigh); MYCHECK(ierr);
+ if(tlow>INT_MAX || thigh>INT_MAX)
+ Msg::Error("Problem too big");
+ *low=tlow;
+ *high=thigh;
}
void LinAlg_GetDoubleInScalar(double *d, gScalar *S)
@@ -403,13 +429,15 @@
void LinAlg_SetScalarInVector(gScalar *S, gVector *V, int i)
{
- ierr = VecSetValues(V->V, 1, &i, &S->s, INSERT_VALUES); MYCHECK(ierr);
+ PetscInt t=i;
+ ierr = VecSetValues(V->V, 1, &t, &S->s, INSERT_VALUES); MYCHECK(ierr);
}
void LinAlg_SetDoubleInVector(double d, gVector *V, int i)
{
PetscScalar tmp = d;
- ierr = VecSetValues(V->V, 1, &i, &tmp, INSERT_VALUES); MYCHECK(ierr);
+ PetscInt t=i;
+ ierr = VecSetValues(V->V, 1, &t, &tmp, INSERT_VALUES); MYCHECK(ierr);
}
void LinAlg_SetComplexInVector(double d1, double d2, gVector *V, int i, int j)
@@ -417,18 +445,22 @@
PetscScalar tmp;
#if PETSC_USE_COMPLEX
tmp = d1 + PETSC_i * d2;
- ierr = VecSetValues(V->V, 1, &i, &tmp, INSERT_VALUES); MYCHECK(ierr);
+ PetscInt ti=i;
+ PetscInt tj=j;
+ ierr = VecSetValues(V->V, 1, &ti, &tmp, INSERT_VALUES); MYCHECK(ierr);
#else
tmp = d1;
- ierr = VecSetValues(V->V, 1, &i, &tmp, INSERT_VALUES); MYCHECK(ierr);
+ ierr = VecSetValues(V->V, 1, &ti, &tmp, INSERT_VALUES); MYCHECK(ierr);
tmp = d2;
- ierr = VecSetValues(V->V, 1, &j, &tmp, INSERT_VALUES); MYCHECK(ierr);
+ ierr = VecSetValues(V->V, 1, &tj, &tmp, INSERT_VALUES); MYCHECK(ierr);
#endif
}
void LinAlg_SetScalarInMatrix(gScalar *S, gMatrix *M, int i, int j)
{
- ierr = MatSetValues(M->M, 1, &i, 1, &j, &S->s, INSERT_VALUES); MYCHECK(ierr);
+ PetscInt ti=i;
+ PetscInt tj=j;
+ ierr = MatSetValues(M->M, 1, &ti, 1, &tj, &S->s, INSERT_VALUES); MYCHECK(ierr);
}
void LinAlg_SetDoubleInMatrix(double d, gMatrix *M, int i, int j)
@@ -455,57 +487,69 @@
void LinAlg_AddScalarInVector(gScalar *S, gVector *V, int i)
{
- ierr = VecSetValues(V->V, 1, &i, &S->s, ADD_VALUES); MYCHECK(ierr);
+ PetscInt ti=i;
+ ierr = VecSetValues(V->V, 1, &ti, &S->s, ADD_VALUES); MYCHECK(ierr);
}
void LinAlg_AddDoubleInVector(double d, gVector *V, int i)
{
PetscScalar tmp = d;
- ierr = VecSetValues(V->V, 1, &i, &tmp, ADD_VALUES); MYCHECK(ierr);
+ PetscInt t=i;
+ ierr = VecSetValues(V->V, 1, &t, &tmp, ADD_VALUES); MYCHECK(ierr);
}
void LinAlg_AddComplexInVector(double d1, double d2, gVector *V, int i, int j)
{
PetscScalar tmp;
+ PetscInt ti=i;
+ PetscInt tj=j;
#if PETSC_USE_COMPLEX
tmp = d1 + PETSC_i * d2;
- ierr = VecSetValues(V->V, 1, &i, &tmp, ADD_VALUES); MYCHECK(ierr);
+ ierr = VecSetValues(V->V, 1, &ti, &tmp, ADD_VALUES); MYCHECK(ierr);
#else
tmp = d1;
- ierr = VecSetValues(V->V, 1, &i, &tmp, ADD_VALUES); MYCHECK(ierr);
+ ierr = VecSetValues(V->V, 1, &ti, &tmp, ADD_VALUES); MYCHECK(ierr);
tmp = d2;
- ierr = VecSetValues(V->V, 1, &j, &tmp, ADD_VALUES); MYCHECK(ierr);
+ ierr = VecSetValues(V->V, 1, &tj, &tmp, ADD_VALUES); MYCHECK(ierr);
#endif
}
void LinAlg_AddScalarInMatrix(gScalar *S, gMatrix *M, int i, int j)
{
- ierr = MatSetValues(M->M, 1, &i, 1, &j, &S->s, ADD_VALUES); MYCHECK(ierr);
+ PetscInt ti=i;
+ PetscInt tj=j;
+ ierr = MatSetValues(M->M, 1, &ti, 1, &tj, &S->s, ADD_VALUES); MYCHECK(ierr);
}
void LinAlg_AddDoubleInMatrix(double d, gMatrix *M, int i, int j)
{
PetscScalar tmp = d;
- ierr = MatSetValues(M->M, 1, &i, 1, &j, &tmp, ADD_VALUES); MYCHECK(ierr);
+ PetscInt ti=i;
+ PetscInt tj=j;
+ ierr = MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, ADD_VALUES); MYCHECK(ierr);
}
void LinAlg_AddComplexInMatrix(double d1, double d2, gMatrix *M, int i, int j, int k, int l)
{
PetscScalar tmp;
+ PetscInt ti=i;
+ PetscInt tj=j;
+ PetscInt tk=k;
+ PetscInt tl=l;
#if PETSC_USE_COMPLEX
tmp = d1 + PETSC_i * d2;
- ierr = MatSetValues(M->M, 1, &i, 1, &j, &tmp, ADD_VALUES); MYCHECK(ierr);
+ ierr = MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, ADD_VALUES); MYCHECK(ierr);
#else
if(d1){
tmp = d1;
- ierr = MatSetValues(M->M, 1, &i, 1, &j, &tmp, ADD_VALUES); MYCHECK(ierr);
- ierr = MatSetValues(M->M, 1, &k, 1, &l, &tmp, ADD_VALUES); MYCHECK(ierr);
+ ierr = MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, ADD_VALUES); MYCHECK(ierr);
+ ierr = MatSetValues(M->M, 1, &tk, 1, &tl, &tmp, ADD_VALUES); MYCHECK(ierr);
}
if(d2){
tmp = -d2;
- ierr = MatSetValues(M->M, 1, &i, 1, &l, &tmp, ADD_VALUES); MYCHECK(ierr);
+ ierr = MatSetValues(M->M, 1, &ti, 1, &tl, &tmp, ADD_VALUES); MYCHECK(ierr);
tmp = d2;
- ierr = MatSetValues(M->M, 1, &k, 1, &j, &tmp, ADD_VALUES); MYCHECK(ierr);
+ ierr = MatSetValues(M->M, 1, &tk, 1, &tj, &tmp, ADD_VALUES); MYCHECK(ierr);
}
#endif
}
@@ -817,7 +861,7 @@
else
Msg::Info("Converged in %d iterations", its);
- for(int i = 0; i < n; i++){
+ for(PetscInt i = 0; i < n; i++){
PetscScalar d;
#if PETSC_USE_COMPLEX
d = solr[i] + PETSC_i * soli[i];
@@ -835,9 +879,9 @@
#endif
-static PetscErrorCode _myKspMonitor(KSP ksp, int it, PetscReal rnorm, void *mctx)
+static PetscErrorCode _myKspMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)
{
- Msg::Info("%3d KSP Residual norm %14.12e", it, rnorm);
+ Msg::Info("%3ld KSP Residual norm %14.12e", (long)it, rnorm);
return 0;
}
@@ -862,9 +906,9 @@
bool view = (!Solver->ksp[kspIndex] && Msg::GetVerbosity() > 0);
if(view && !RankCpu){
- int i, j;
+ PetscInt i, j;
ierr = MatGetSize(A->M, &i, &j); MYCHECK(ierr);
- Msg::Info("N: %d", i);
+ Msg::Info("N: %ld", (long)i);
}
if(kspIndex != 0){
@@ -908,7 +952,7 @@
}
if(!RankCpu){
- int its;
+ PetscInt its;
ierr = KSPGetIterationNumber(Solver->ksp[kspIndex], &its); MYCHECK(ierr);
Msg::Info("%d iterations", its);
}