linalg_sparskit.c

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

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

  GetDP_Begin("LinAlg_GetDoubleInMatrix");

  get_element_in_matrix(&M->M, i, j, d);

  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_GetScalarInMatrix' not yet implemented");  

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_GetColumnInMatrix");

  get_column_in_matrix(&M->M, col, V1->V);

  GetDP_End ;
}

void LinAlg_GetMatrixContext(gMatrix *A, void **ctx){

  GetDP_Begin("LinAlg_GetMatrixContext");

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

  GetDP_End;
}

/* Set */

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

  GetDP_Begin("LinAlg_SetScalar");

  S->d = d[0] ;

  GetDP_End ;
}

void LinAlg_SetVector(gVector *V, double *v){
  int i;

  GetDP_Begin("LinAlg_SetScalar");

  for(i=0; i<V->N; i++) V->V[i] = *v ;
  
  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_SetScalarInVector");

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

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_SetDoubleInVector");

  V->V[i] = d ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_SetComplexInVector");

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

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_SetScalarInMatrix");

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

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_SetDoubleInMatrix");

  Msg(GERROR, "'LinAlg_SetScalarInMatrix' 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_SetScalarInMatrix' not yet implemented");  

  GetDP_End ;
}

/* Add */

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

  GetDP_Begin("LinAlg_AddScalarScalar");

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

  GetDP_End ;
}

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

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

  for (i=0 ; i<V->N ; i++) if (DummyDof[i] == 1) V->V[i] = 0 ;
  GetDP_End ;
}

void LinAlg_AddScalarInVector(gScalar *S, gVector *V, int i){
  int * DummyDof;
  GetDP_Begin("LinAlg_AddScalarInVector");

  if ((DummyDof = Current.DofData->DummyDof)) 
    if (DummyDof[i] == 1) GetDP_End ;

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

  GetDP_End ;
}

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

  int * DummyDof;
  GetDP_Begin("LinAlg_AddDoubleInVector");

  if ((DummyDof = Current.DofData->DummyDof))
    if (DummyDof[i] == 1)
      GetDP_End ;

  V->V[i] += d ;

  GetDP_End ;
}

void LinAlg_AddComplexInVector(double d1, double d2, gVector *V, int i, int j){
  int * DummyDof, iok,jok;
  GetDP_Begin("LinAlg_AddComplexInVector");

  iok=jok=1;

  if ((DummyDof = Current.DofData->DummyDof)) { 
    if (DummyDof[i] == 1) iok=0;
    if (DummyDof[j] == 1) jok=0;
  }

  if (iok)
    V->V[i] += d1 ;

  if (jok)
    V->V[j] += d2 ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_AddScalarInMatrix");

  add_matrix_double(&M->M, i+1, j+1, S->d) ;

  GetDP_End ;
}

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

  int * DummyDof;

  GetDP_Begin("LinAlg_AddDoubleInMatrix");

  if ((DummyDof = Current.DofData->DummyDof))
    if ( (DummyDof[i] == 1 || DummyDof[j] == 1) && (i != j) ) 
    GetDP_End ;

  add_matrix_double(&M->M, i+1, j+1, d) ;

  GetDP_End ;
}

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

  GetDP_Begin("LinAlg_AddComplexInMatrix");

  if(d1){
    add_matrix_double(&M->M, i+1, j+1, d1) ;
    add_matrix_double(&M->M, k+1, l+1, d1) ;
  }
  if(d2){
    add_matrix_double(&M->M, i+1, l+1, -d2) ;
    add_matrix_double(&M->M, k+1, j+1, d2) ;
  }

  GetDP_End ;
}

void LinAlg_AddVectorVector(gVector *V1, gVector *V2, gVector *V3){

  GetDP_Begin("LinAlg_AddvectorVector");

  if(V3 == V1)
    add_vector_vector(V1->N, V1->V, V2->V) ;
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_AddVectorVector'");

  GetDP_End ;
}

void LinAlg_AddVectorProdVectorDouble(gVector *V1, gVector *V2, double d, gVector *V3){

  GetDP_Begin("LinAlg_AddVectorProdVectorDouble");

  if(V3 == V1)
    add_vector_prod_vector_double(V1->N, V1->V, V2->V, d) ;
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_AddVectorProdVectorDouble'");

  GetDP_End ;
}

void LinAlg_AddMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3){

  GetDP_Begin("LinAlg_AddMatrixMatrix");

  if(M3 == M1)
    add_matrix_matrix(&M1->M, &M2->M) ;
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_AddMatrixMatrix'");

  GetDP_End ;
}

void LinAlg_AddMatrixProdMatrixDouble(gMatrix *M1, gMatrix *M2, double d, gMatrix *M3){

  GetDP_Begin("LinAlg_AddMatrixProdMatrixDouble");

  if(M3 == M1)
    add_matrix_prod_matrix_double(&M1->M, &M2->M, d) ;
  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->d = S1->d - S2->d ;

  GetDP_End ;
}

void LinAlg_SubVectorVector(gVector *V1, gVector *V2, gVector *V3){

  GetDP_Begin("LinAlg_SubVectorVector");

  if(V3 == V1) 
    sub_vector_vector_1(V1->N, V1->V, V2->V) ;
  else if (V3 == V2) 
    sub_vector_vector_2(V1->N, V1->V, V2->V) ;
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_SubVectorVector'");  

  GetDP_End ;
}

void LinAlg_SubMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3){

  GetDP_Begin("LinAlg_SubMatrixMatrix");

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

  GetDP_End ;
}

/* Prod */

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

  GetDP_Begin("LinAlg_ProdScalarScalar");

  S3->d = S1->d * S2->d ;

  GetDP_End ;
}

void LinAlg_ProdScalarDouble(gScalar *S1, double d, gScalar *S2){

  GetDP_Begin("LinAlg_ProdScalarDouble");

  S2->d = S1->d * d ;

  GetDP_End ;
}

void LinAlg_ProdScalarComplex(gScalar *S, double d1, double d2, double *d3, double *d4){

  GetDP_Begin("LinAlg_ProdScalarComplex");

  *d3 = S->d * d1 ;
  *d4 = S->d * d2 ;

  GetDP_End ;
}

void LinAlg_ProdVectorScalar(gVector *V1, gScalar *S, gVector *V2){

  GetDP_Begin("LinAlg_ProdvectorScalar");

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

  GetDP_End ;
}

void LinAlg_ProdVectorDouble(gVector *V1, double d, gVector *V2){
  
  GetDP_Begin("LinAlg_ProdVectorDouble");

  if(V2 == V1)
    prod_vector_double(V1->N, V1->V, d);
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdVectorDouble'");

  GetDP_End ;
}

void LinAlg_ProdVectorComplex(gVector *V1, double d1, double d2, gVector *V2){

  GetDP_Begin("LinAlg_ProdVectorComplex");

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

  GetDP_End ;
}

void LinAlg_ProdVectorVector(gVector *V1, gVector *V2, double *d){

  GetDP_Begin("LinAlg_ProdVectorVector");

  prodsc_vector_vector (V1->N, V1->V, V2->V, d) ;

  GetDP_End ;
}

void LinAlg_ProdMatrixVector(gMatrix *M, gVector *V1, gVector *V2){

  GetDP_Begin("LinAlg_ProdMatrixVector");

  if(V2 == V1)
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdMatrixVector'");
  else
    prod_matrix_vector(&M->M, V1->V, V2->V);

  GetDP_End ;
}

void LinAlg_ProdMatrixScalar(gMatrix *M1, gScalar *S, gMatrix *M2){

  GetDP_Begin("LinAlg_ProdMatrixScalar");

  if(M2 == M1)
    prod_matrix_double (&M1->M, S->d);
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdMatrixScalar'");

  GetDP_End ;
}

void LinAlg_ProdMatrixDouble(gMatrix *M1, double d, gMatrix *M2){

  GetDP_Begin("LinAlg_ProdMatrixDouble");

  if(M2 == M1)
    prod_matrix_double (&M1->M, d);
  else
    Msg(GERROR, "Wrong arguments in 'LinAlg_ProdMatrixDouble'");

  GetDP_End ;
}

void LinAlg_ProdMatrixComplex(gMatrix *M1, double d1, double d2, gMatrix *M2){

  GetDP_Begin("LinAlg_ProdMatrixComplex");

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

  GetDP_End ;
}

/* Div */

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

  GetDP_Begin("LinAlg_DivScalarScalar");

  S3->d = S1->d / S2->d ;

  GetDP_End ;
}

void LinAlg_DivScalarDouble(gScalar *S1, double d, gScalar *S2){

  GetDP_Begin("LinAlg_DivScalarDouble");

  S2->d = S1->d / d ;

  GetDP_End ;
}

/* Norm */

void LinAlg_VectorNorm2(gVector *V1, double *norm){

  GetDP_Begin("LinAlg_VectorNorm2");

  norm2_vector(V1->N, V1->V, norm);
  
  GetDP_End ;
}

void LinAlg_VectorNormInf(gVector *V1, double *norm){

  GetDP_Begin("LinAlg_VectorNormInf");
  
  norminf_vector(V1->N, V1->V, norm);

  GetDP_End;
}

/* Assemble */

void LinAlg_AssembleMatrix(gMatrix *M){
}

void LinAlg_AssembleVector(gVector *V){
}

/* FMM */

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

  GetDP_Begin("LinAlg_MatVectorProdVector");

  FMM_MatVectorProd(V1->V, V2->V) ;

  GetDP_End ;
}

/* Solve */

void LinAlg_Solve(gMatrix *A, gVector *B, gSolver *Solver, gVector *X){

  GetDP_Begin("LinAlg_Solve");

  solve_matrix(&A->M, &Solver->Params, B->V, X->V);

  GetDP_End ;
}

void LinAlg_SolveAgain(gMatrix *A, gVector *B, gSolver *Solver, gVector *X){
  int tmp;

  GetDP_Begin("LinAlg_SolveAgain");

  tmp = Solver->Params.Re_Use_LU;
  Solver->Params.Re_Use_LU = 1;
  solve_matrix(&A->M, &Solver->Params, B->V, X->V);
  Solver->Params.Re_Use_LU = tmp;

  GetDP_End ;
}

#endif