linalg_petsc.c

cfd求解器使用与gmsh网格的求解

/* Write */

void LinAlg_WriteScalar(FILE *file, gScalar *S){

  GetDP_Begin("LinAlg_WriteScalar");

  Msg(GERROR, "'LinAlg_WriteScalar' not yet implemented");

  GetDP_End;
}

void LinAlg_WriteVector(FILE *file, gVector *V){
  PetscScalar *tmp;
  int n;

  GetDP_Begin("LinAlg_WriteVector");

  ierr = VecGetLocalSize(V->V, &n); MYCHECK(ierr);
  ierr = VecGetArray(V->V, &tmp); MYCHECK(ierr);
  safe_fwrite(tmp, sizeof(PetscScalar), n, file);
  fprintf(file, "\n");
  ierr = VecRestoreArray(V->V, &tmp); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_WriteMatrix(FILE *file, gMatrix *M){

  GetDP_Begin("LinAlg_WriteMatrix");

  Msg(GERROR, "'LinAlg_WriteMatrix' not yet implemented");  

  GetDP_End;
}

/* Get */

void LinAlg_GetVectorSize(gVector *V, int *i){

  GetDP_Begin("LinAlg_GetVectorSize");

  ierr = VecGetSize(V->V, i); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_GetLocalVectorRange(gVector *V, int *low, int *high){

  GetDP_Begin("LinAlg_GetLocalVectorRange");

  ierr = VecGetOwnershipRange(V->V, low, high); MYCHECK(ierr);

  GetDP_End;

}

void LinAlg_GetMatrixSize(gMatrix *M, int *i, int *j){

  GetDP_Begin("LinAlg_GetMatrixSize");

  ierr = MatGetSize(M->M, i, j); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_GetLocalMatrixRange(gMatrix *M, int *low, int *high){

  GetDP_Begin("LinAlg_GetLocalMatrixRange");

  ierr = MatGetOwnershipRange(M->M, low, high); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_GetDoubleInScalar(double *d, gScalar *S){

  GetDP_Begin("LinAlg_GetDoubleInScalar");

#if PETSC_USE_COMPLEX
  *d = real(S->s);
#else
  *d = S->s;
#endif

  GetDP_End;
}

void LinAlg_GetComplexInScalar(double *d1, double *d2, gScalar *S){

  GetDP_Begin("LinAlg_GetComplexInScalar");

#if PETSC_USE_COMPLEX
  *d1 = real(S->s);
  *d2 = imag(S->s);
#else
  Msg(GERROR, "'LinAlg_GetComplexInScalar' not available with this Solver");  
#endif

  GetDP_End;
}

void LinAlg_GetScalarInVector(gScalar *S, gVector *V, int i){
  PetscScalar *tmp;

  GetDP_Begin("LinAlg_GetScalarInVector");

  ierr = VecGetArray(V->V, &tmp); MYCHECK(ierr);
  S->s = tmp[i];
  ierr = VecRestoreArray(V->V, &tmp); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_GetDoubleInVector(double *d, gVector *V, int i){
  PetscScalar *tmp;

  GetDP_Begin("LinAlg_GetDoubleInVector");

  ierr = VecGetArray(V->V, &tmp); MYCHECK(ierr);
#if PETSC_USE_COMPLEX
  *d = real(tmp[i]);
#else
  *d = tmp[i];
#endif
  ierr = VecRestoreArray(V->V, &tmp); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_GetAbsDoubleInVector(double *d, gVector *V, int i){
  PetscScalar *tmp;

  GetDP_Begin("LinAlg_GetAbsDoubleInVector");

  ierr = VecGetArray(V->V, &tmp); MYCHECK(ierr);
#if PETSC_USE_COMPLEX
  *d = fabs(real(tmp[i]));
#else
  *d = fabs(tmp[i]);
#endif
  ierr = VecRestoreArray(V->V, &tmp); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_GetComplexInVector(double *d1, double *d2, gVector *V, int i, int j){
  PetscScalar *tmp;

  GetDP_Begin("LinAlg_GetComplexInVector");

  ierr = VecGetArray(V->V, &tmp); MYCHECK(ierr);
#if PETSC_USE_COMPLEX
  *d1 = real(tmp[i]);
  *d2 = imag(tmp[i]);
#else
  *d1 = (double)tmp[i];
  *d2 = (double)tmp[j];
#endif
  ierr = VecRestoreArray(V->V, &tmp); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_GetScalarInMatrix(gScalar *S, gMatrix *M, int i, int j){

  GetDP_Begin("LinAlg_GetScalarInMatrix");

  Msg(GERROR, "'LinAlg_GetScalarInMatrix' not yet implemented");  

  GetDP_End;
}

void LinAlg_GetDoubleInMatrix(double *d, gMatrix *M, int i, int j){

  GetDP_Begin("LinAlg_GetDoubleInMatrix");

  Msg(GERROR, "'LinAlg_GetDoubleInMatrix' not yet implemented");  

  GetDP_End;
}

void LinAlg_GetComplexInMatrix(double *d1, double *d2, gMatrix *M, int i, int j, int k, int l){

  GetDP_Begin("LinAlg_GetComplexInMatrix");

  Msg(GERROR, "'LinAlg_GetComplexInMatrix' not yet implemented");  

  GetDP_End;
}

void LinAlg_GetColumnInMatrix(gMatrix *M, int col, gVector *V1){

  GetDP_Begin("LinAlg_GetColumnInMatrix");

  Msg(GERROR, "'LinAlg_GetColumnInMatrix' not yet implemented");  

  GetDP_End;
}

void LinAlg_GetMatrixContext(gMatrix *A, void **myCtx){
  petscCtx *ctx;

  GetDP_Begin("LinAlg_GetMatrixContext");

  MatShellGetContext(A->M,(void**)&ctx);
  *myCtx = ctx->ctx;
  
  GetDP_End;
}

/* Set */

void LinAlg_SetScalar(gScalar *S, double *d){

  GetDP_Begin("LinAlg_SetScalar");

#if PETSC_USE_COMPLEX
  S->s = d[0] + (PETSC_i * d[1]);
#else
  S->s = d[0];
#endif

  GetDP_End;
}

void LinAlg_SetVector(gVector *V, double *v){
  PetscScalar tmp = *v;

  GetDP_Begin("LinAlg_SetVector");

  ierr = VecSet(V->V, tmp); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_SetScalarInVector(gScalar *S, gVector *V, int i){

  GetDP_Begin("LinAlg_SetScalarInVector");

  ierr = VecSetValues(V->V, 1, &i, &S->s, INSERT_VALUES); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_SetDoubleInVector(double d, gVector *V, int i){
  PetscScalar tmp;

  GetDP_Begin("LinAlg_SetDoubleInVector");

  tmp = d;
  ierr = VecSetValues(V->V, 1, &i, &tmp, INSERT_VALUES); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_SetComplexInVector(double d1, double d2, gVector *V, int i, int j){
  PetscScalar tmp;

  GetDP_Begin("LinAlg_SetComplexInVector");

#if PETSC_USE_COMPLEX
  tmp = d1 + PETSC_i * d2;
  ierr = VecSetValues(V->V, 1, &i, &tmp, INSERT_VALUES); MYCHECK(ierr);
#else
  tmp = d1;
  ierr = VecSetValues(V->V, 1, &i, &tmp, INSERT_VALUES); MYCHECK(ierr);
  tmp = d2;
  ierr = VecSetValues(V->V, 1, &j, &tmp, INSERT_VALUES); MYCHECK(ierr);
#endif

  GetDP_End;
}

void LinAlg_SetScalarInMatrix(gScalar *S, gMatrix *M, int i, int j){

  GetDP_Begin("LinAlg_SetScalarInMatrix");

  ierr = MatSetValues(M->M, 1, &i, 1, &j, &S->s, INSERT_VALUES); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_SetDoubleInMatrix(double d, gMatrix *M, int i, int j){

  GetDP_Begin("LinAlg_SetDoubleInMatrix");

  Msg(GERROR, "'LinAlg_SetDoubleInMatrix' not yet implemented");  

  GetDP_End;
}

void LinAlg_SetComplexInMatrix(double d1, double d2, gMatrix *M, int i, int j, int k, int l){

  GetDP_Begin("LinAlg_SetComplexInMatrix");

  Msg(GERROR, "'LinAlg_SetComplexInMatrix' not yet implemented");  

  GetDP_End;
}

/* Add */

void LinAlg_AddScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3){

  GetDP_Begin("LinAlg_AddScalarScalar");

  S3->s = S1->s + S2->s;

  GetDP_End;
}

void LinAlg_DummyVector(gVector *V){
  int * DummyDof;
  
  GetDP_Begin("LinAlg_DummyVector");

  DummyDof = Current.DofData->DummyDof;
  if(DummyDof == NULL) GetDP_End;

  Msg(GERROR, "'LinAlg_DummyVector' not yet implemented");  

  GetDP_End;
}

void LinAlg_AddScalarInVector(gScalar *S, gVector *V, int i){

  GetDP_Begin("LinAlg_AddScalarInVector");

  ierr = VecSetValues(V->V, 1, &i, &S->s, ADD_VALUES); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_AddDoubleInVector(double d, gVector *V, int i){
  PetscScalar tmp;

  GetDP_Begin("LinAlg_AddDoubleInVector");

  tmp = d;
  ierr = VecSetValues(V->V, 1, &i, &tmp, ADD_VALUES); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_AddComplexInVector(double d1, double d2, gVector *V, int i, int j){
  PetscScalar tmp;

  GetDP_Begin("LinAlg_AddComplexInVector");

#if PETSC_USE_COMPLEX
  tmp = d1 + PETSC_i * d2;
  ierr = VecSetValues(V->V, 1, &i, &tmp, ADD_VALUES); MYCHECK(ierr);
#else
  tmp = d1;
  ierr = VecSetValues(V->V, 1, &i, &tmp, ADD_VALUES); MYCHECK(ierr);
  tmp = d2;
  ierr = VecSetValues(V->V, 1, &j, &tmp, ADD_VALUES); MYCHECK(ierr);
#endif

  GetDP_End;
}

void LinAlg_AddScalarInMatrix(gScalar *S, gMatrix *M, int i, int j){

  GetDP_Begin("LinAlg_AddScalarInMatrix");

  ierr = MatSetValues(M->M, 1, &i, 1, &j, &S->s, ADD_VALUES); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_AddDoubleInMatrix(double d, gMatrix *M, int i, int j){
  PetscScalar tmp;

  GetDP_Begin("LinAlg_AddDoubleInMatrix");

  tmp = d;
  ierr = MatSetValues(M->M, 1, &i, 1, &j, &tmp, ADD_VALUES); MYCHECK(ierr);

  GetDP_End;
}

void LinAlg_AddComplexInMatrix(double d1, double d2, gMatrix *M, int i, int j, int k, int l){
  PetscScalar tmp;

  GetDP_Begin("LinAlg_AddComplexInMatrix");

#if PETSC_USE_COMPLEX
  tmp = d1 + PETSC_i * d2;
  ierr = MatSetValues(M->M, 1, &i, 1, &j, &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);
  }
  if(d2){
    tmp = -d2;
    ierr = MatSetValues(M->M, 1, &i, 1, &l, &tmp, ADD_VALUES); MYCHECK(ierr);
    tmp = d2;
    ierr = MatSetValues(M->M, 1, &k, 1, &j, &tmp, ADD_VALUES); MYCHECK(ierr);
  }
#endif

  GetDP_End;
}

void LinAlg_AddVectorVector(gVector *V1, gVector *V2, gVector *V3){
  PetscScalar tmp = 1.0;

  GetDP_Begin("LinAlg_AddVectorVector");

  if(V3 == V1){
    ierr = VecAXPY(V1->V, tmp, V2->V); MYCHECK(ierr);
  }
  else if(V3 == V2){
    ierr = VecAXPY(V2->V, tmp, V1->V); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_AddVectorVector'");  

  GetDP_End;
}

void LinAlg_AddVectorProdVectorDouble(gVector *V1, gVector *V2, double d, gVector *V3){
  PetscScalar tmp = d;

  GetDP_Begin("LinAlg_AddvectorProdVectorDouble");

  if(V3 == V1){
    ierr = VecAXPY(V1->V, tmp, V2->V); MYCHECK(ierr);
  }
  else if(V3 == V2){
    ierr = VecAYPX(V2->V, tmp, V1->V); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_AddVectorProdVectorDouble'");  

  GetDP_End;
}

void LinAlg_AddMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3){
  PetscScalar tmp = 1.0;

  GetDP_Begin("LinAlg_AddMatrixMatrix");

  if(M3 == M1){
    ierr = MatAXPY(M1->M, tmp, M2->M, DIFFERENT_NONZERO_PATTERN); MYCHECK(ierr);
  }
  else if(M3 == M2){
    ierr = MatAXPY(M2->M, tmp, M1->M, DIFFERENT_NONZERO_PATTERN); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_AddMatrixMatrix'");  

  GetDP_End;
}

void LinAlg_AddMatrixProdMatrixDouble(gMatrix *M1, gMatrix *M2, double d, gMatrix *M3){
  PetscScalar tmp = d;

  GetDP_Begin("LinAlg_AddMatrixProdMatrixDouble");

  if(M3 == M1){
    ierr = MatAXPY(M1->M, tmp, M2->M, DIFFERENT_NONZERO_PATTERN); MYCHECK(ierr);
  }
  else if(M3 == M2){
    ierr = MatAYPX(M2->M, tmp, M1->M); MYCHECK(ierr);
  }
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_AddMatrixProdMatrixDouble'");

  GetDP_End;
}

/* Sub */

void LinAlg_SubScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3){

  GetDP_Begin("LinAlg_SubScalarScalar");

  S3->s = S1->s - S2->s;

  GetDP_End;
}

void LinAlg_SubVectorVector(gVector *V1, gVector *V2, gVector *V3){
  PetscScalar tmp = -1.0;

  GetDP_Begin("LinAlg_SubVectorVector");

  if(V3 == V1){
    ierr = VecAXPY(V1->V, tmp, V2->V); MYCHECK(ierr);
  }