linalg_sparskit.c

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

#define RCSID "$Id: LinAlg_SPARSKIT.c,v 1.35 2006/02/26 00:42:53 geuzaine Exp $"
/*
 * Copyright (C) 1997-2006 P. Dular, C. Geuzaine
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA.
 *
 * Please report all bugs and problems to <[e-mail]>.
 *
 * Contributor(s):
 *   Ruth Sabariego
 */

#if defined(HAVE_SPARSKIT)

/* This is the interface library for the Sparskit solvers */

#include "GetDP.h"
#include "Magic.h"
#include "LinAlg.h"
#include "SafeIO.h"
#include "F_FMMOperations.h"

#include "Data_Passive.h"
#include "CurrentData.h"

extern char  Name_Path[MAX_FILE_NAME_LENGTH] ;
static char *Name_SolverFile = NULL, *Name_DefaultSolverFile = "solver.par" ;
static char *SolverOptions[100];

/* Init */

void LinAlg_Initialize(int* argc, char*** argv, int *NbrCpu, int *RankCpu){

  GetDP_Begin("LinAlg_Initialize");

  *NbrCpu = 1 ;
  *RankCpu = 0 ;

  GetDP_End ;
}

void LinAlg_InitializeSolver(int* sargc, char*** sargv, int *NbrCpu, int *RankCpu){
  int i=1, argc, iopt=0;
  char **argv;
 
  GetDP_Begin("LinAlg_InitializeSolver");

  *NbrCpu = 1 ;
  *RankCpu = 0 ;
  argc = *sargc ;
  argv = *sargv ;
  SolverOptions[0] = NULL;

  while (i < argc) {
    if (argv[i][0] == '-') {
      if (!strcmp(argv[i]+1, "solver") || !strcmp(argv[i]+1, "s")) {
	i++ ;
	if (i<argc && argv[i][0]!='-')
	  Name_SolverFile = argv[i++] ;
	else
	  Msg(GERROR, "Missing file name");
      }
      else{
	i++ ;
	if (i<argc && argv[i][0]!='-') { 
	  SolverOptions[iopt++] = argv[i-1]+1;
	  SolverOptions[iopt++] = argv[i];
	  SolverOptions[iopt] = NULL;
	  i++ ; 
	}
	else {
	  Msg(GERROR, "Missing number");
	}
      }
    }
    else
      Msg(SPARSKIT, "Unknown option: '%s'\n", argv[i++]) ; 
  }

  GetDP_End ;
}

/* Finalize */

void LinAlg_Finalize(void){
}

void LinAlg_FinalizeSolver(void){
}

/* Barrier */

void LinAlg_Barrier(void){
}

/* Sequential */

void LinAlg_SequentialBegin(void){
}

void LinAlg_SequentialEnd(void){
}

/* Create */

void LinAlg_CreateSolver(gSolver *Solver, char * SolverDataFileName){
  int  i;
  char FileName[MAX_FILE_NAME_LENGTH];

  GetDP_Begin("LinAlg_CreateSolver");
  
  strcpy(FileName, Name_Path);

  if(SolverDataFileName){
    /* name in .pro file */
    if(SolverDataFileName[0] == '/' || SolverDataFileName[0] == '\\'){
      /* -> absolute if it starts with / or \ */
      strcpy(FileName, SolverDataFileName);
    }
    else{
      /* -> relative otherwise */
      strcat(FileName, SolverDataFileName);
    }
  }
  else if (Name_SolverFile){
    /* name on command line -> always absolute */
    strcpy(FileName, Name_SolverFile);
  }
  else{
    /* default file name -> always relative */
    strcat(FileName, Name_DefaultSolverFile);
  }
  
  Msg(SPARSKIT, "Loading parameter file '%s'\n", FileName);

  init_solver(&Solver->Params, FileName) ;

  i = 0;
  while(SolverOptions[i] && SolverOptions[i+1]){
    init_solver_option(&Solver->Params, SolverOptions[i], SolverOptions[i+1]) ;
    i+=2;
  }

  GetDP_End ;
}

void LinAlg_CreateVector(gVector *V, gSolver *Solver, int n, int NbrPart, int *Part){

  GetDP_Begin("LinAlg_CreateVector");

  init_vector(n, &V->V) ;
  V->N = n ;

  GetDP_End ;
}

void LinAlg_CreateMatrix(gMatrix *M, gSolver *Solver, int n, int m, 
			 int NbrPart, int *Part, int *Nnz){

  GetDP_Begin("LinAlg_CreateMatrix");

  init_matrix(n, &M->M, &Solver->Params) ;

  GetDP_End ;
}

void LinAlg_CreateMatrixShell(gMatrix *A, gSolver *Solver, int n, int m, void *ctx, 
			      void (*myMult)(gMatrix *A, gVector *x, gVector *y)){
  GetDP_Begin("LinAlg_CreateMatrixShell");

  Msg(GERROR, "Matrix-free operations not implemented with Sparskit");

  GetDP_End;
}

/* Destroy */

void LinAlg_DestroySolver(gSolver *Solver){

  GetDP_Begin("LinAlg_DestroySolver");

  Msg(WARNING, "'LinAlg_DestroySolver' not yet implemented");

  GetDP_End ;
}

void LinAlg_DestroyVector(gVector *V){

  GetDP_Begin("LinAlg_DestroyVector");

  Free(V->V);

  GetDP_End ;
}

void LinAlg_DestroyMatrix(gMatrix *M){

  GetDP_Begin("LinAlg_DestroyMatrix");

  Msg(WARNING, "'LinAlg_DestroyMatrix' not yet implemented");  

  GetDP_End ;
}

/* Copy */

void LinAlg_CopyScalar(gScalar *S1, gScalar *S2){

  GetDP_Begin("LinAlg_CopyScalar");

  S2->d = S1->d ;

  GetDP_End ;
}

void LinAlg_CopyVector(gVector *V1, gVector *V2){

  GetDP_Begin("LinAlg_CopyVector");

  memcpy(V2->V, V1->V, V1->N*sizeof(double)) ;

  GetDP_End ;
}

void LinAlg_CopyMatrix(gMatrix *M1, gMatrix *M2){

  GetDP_Begin("LinAlg_CopyMatrix");

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

  GetDP_End ;
}

/* Zero */

void LinAlg_ZeroScalar(gScalar *S){

  GetDP_Begin("LinAlg_ZeroScalar");

  S->d = 0. ;

  GetDP_End ;
}

void LinAlg_ZeroVector(gVector *V){

  GetDP_Begin("LinAlg_ZeroCopyVector");

  zero_vector(V->N, V->V) ;

  GetDP_End ;
}

void LinAlg_ZeroMatrix(gMatrix *M){
  int i ;

  GetDP_Begin("LinAlg_ZeroMatrix");

  zero_matrix(&M->M) ;
  /* la routine de produit matrice vecteur est buggee s'il
     existe des lignes sans aucun element dans la matrice... */
  for(i=0 ; i<M->M.N ; i++) add_matrix_double(&M->M, i+1, i+1, 0.0) ;

  GetDP_End ;
}

/* Scan */

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

  GetDP_Begin("LinAlg_ScanScalar");

  fscanf(file, "%lf", &S->d) ;

  GetDP_End ;
}

void LinAlg_ScanVector(FILE *file, gVector *V) {
  int i ;

  GetDP_Begin("LinAlg_ScanVector");

  for(i=0 ; i<V->N ; i++) fscanf(file, "%lf", &V->V[i]) ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_ScanMatrix");

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

  GetDP_End ;
}

/* Read */

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

  GetDP_Begin("LinAlg_ReadScalar");

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

  GetDP_End ;
}

void LinAlg_ReadVector(FILE *file, gVector *V) {

  GetDP_Begin("LinAlg_ReadVector");

  fread(V->V, sizeof(double), V->N, file);

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_ReadMatrix");

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

  GetDP_End ;
}

/* Print */

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

  GetDP_Begin("LinAlg_PrintScalar");

  fprintf(file, "%.16g", S->d) ;

  GetDP_End ;
}

void LinAlg_PrintVector(FILE *file, gVector *V){

  GetDP_Begin("LinAlg_PrintVector");

  formatted_write_vector(file, V->N, V->V, KUL) ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_PrintMatrix");

  formatted_write_matrix(file, &M->M, KUL) ;

  GetDP_End ;
}

/* 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){

  GetDP_Begin("LinAlg_WriteVector");

  safe_fwrite(V->V, sizeof(double), V->N, file);
  fprintf(file, "\n");

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_WriteMatrix");

  binary_write_matrix(&M->M, "A", ".mat") ;

  GetDP_End ;
}

/* Get */

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

  GetDP_Begin("LinAlg_GetVectorSize");

  *i = V->N ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetLocalVectorRange");

  *low = 0 ;
  *high = V->N ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetMatrixSize");

  *i = *j = M->M.N ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetLocalVectorRange");

  *low = 0 ;
  *high = M->M.N ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetDoubleInScalar");

  *d = S->d ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetComplexInScalar");

  Msg(GERROR, "'LinAlg_GetComplexInScalar' not available with this Solver");

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetScalarInVector");

  S->d = V->V[i] ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetDoubleInVector");

  *d = V->V[i] ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetAbsDoubleInVector");

  *d = fabs(V->V[i]) ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetComplexInVector");

  *d1 = V->V[i] ;
  *d2 = V->V[j] ;

  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 ;
}