is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.
| LinAlg_SPARSKIT.cpp (getdp-3.4.0-source.tgz) | : | LinAlg_SPARSKIT.cpp (getdp-3.5.0-source.tgz) |
|---|
|
| // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege | | // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege |
| // | | // |
| // See the LICENSE.txt file for license information. Please report all | | // See the LICENSE.txt file for license information. Please report all |
| // issues on https://gitlab.onelab.info/getdp/getdp/issues. | | // issues on https://gitlab.onelab.info/getdp/getdp/issues. |
| // | | // |
| // Contributor(s): | | // Contributor(s): |
| // Ruth Sabariego | | // Ruth Sabariego |
| // | | // |
| | |
| #include <math.h> | | #include <math.h> |
| #include <string.h> | | #include <string.h> |
| #include "GetDPConfig.h" | | #include "GetDPConfig.h" |
| #include "LinAlg.h" | | #include "LinAlg.h" |
| #include "ProData.h" | | #include "ProData.h" |
| #include "DofData.h" | | #include "DofData.h" |
| #include "MallocUtils.h" | | #include "MallocUtils.h" |
| #include "Message.h" | | #include "Message.h" |
| | |
|
| extern struct CurrentData Current ; | | extern struct CurrentData Current; |
| | |
| extern char *Name_Path; | | extern char *Name_Path; |
| | |
| #if defined(HAVE_SPARSKIT) | | #if defined(HAVE_SPARSKIT) |
| | |
|
| static const char *Name_SolverFile = NULL, *Name_DefaultSolverFile = "solver.par | | static const char *Name_SolverFile = NULL, |
| " ; | | *Name_DefaultSolverFile = "solver.par"; |
| static char *SolverOptions[100]; | | static char *SolverOptions[100]; |
| | |
|
| void LinAlg_InitializeSolver(int* sargc, char*** sargv) | | void LinAlg_InitializeSolver(int *sargc, char ***sargv) |
| { | | { |
|
| int i=1, argc, iopt=0; | | int i = 1, argc, iopt = 0; |
| char **argv; | | char **argv; |
| | |
|
| argc = *sargc ; | | argc = *sargc; |
| argv = *sargv ; | | argv = *sargv; |
| SolverOptions[0] = NULL; | | SolverOptions[0] = NULL; |
| | |
|
| while (i < argc) { | | while(i < argc) { |
| if (argv[i][0] == '-') { | | if(argv[i][0] == '-') { |
| if (!strcmp(argv[i]+1, "solver") || !strcmp(argv[i]+1, "s")) { | | if(!strcmp(argv[i] + 1, "solver") || !strcmp(argv[i] + 1, "s")) { |
| i++ ; | | i++; |
| if (i<argc && argv[i][0]!='-') | | if(i < argc && argv[i][0] != '-') |
| Name_SolverFile = argv[i++] ; | | Name_SolverFile = argv[i++]; |
| else | | else |
| Message::Error("Missing file name"); | | Message::Error("Missing file name"); |
| } | | } |
|
| else if (!strcmp(argv[i]+1, "slepc")) { | | else if(!strcmp(argv[i] + 1, "slepc")) { |
| i++; | | i++; |
| } | | } |
|
| else{ | | else { |
| i++ ; | | i++; |
| if (i<argc && argv[i][0]!='-') { | | if(i < argc && argv[i][0] != '-') { |
| SolverOptions[iopt++] = argv[i-1]+1; | | SolverOptions[iopt++] = argv[i - 1] + 1; |
| SolverOptions[iopt++] = argv[i]; | | SolverOptions[iopt++] = argv[i]; |
| SolverOptions[iopt] = NULL; | | SolverOptions[iopt] = NULL; |
| i++ ; | | i++; |
| } | | } |
| else { | | else { |
| Message::Error("Missing number"); | | Message::Error("Missing number"); |
| } | | } |
| } | | } |
| } | | } |
| else | | else |
|
| Message::Info("Unknown option: '%s'", argv[i++]) ; | | Message::Info("Unknown option: '%s'", argv[i++]); |
| } | | } |
| } | | } |
| | |
|
| void LinAlg_FinalizeSolver() | | void LinAlg_FinalizeSolver() {} |
| { | | |
| } | | |
| | |
|
| void LinAlg_SetCommSelf() | | void LinAlg_SetCommSelf() {} |
| { | | |
| } | | |
| | |
|
| void LinAlg_SetCommWorld() | | void LinAlg_SetCommWorld() {} |
| { | | |
| } | | |
| | |
| void LinAlg_CreateSolver(gSolver *Solver, const char *SolverDataFileName) | | void LinAlg_CreateSolver(gSolver *Solver, const char *SolverDataFileName) |
| { | | { |
|
| int i; | | int i; |
| char FileName[256]; | | char FileName[256]; |
| | |
| strcpy(FileName, Name_Path); | | strcpy(FileName, Name_Path); |
| | |
|
| if(SolverDataFileName){ | | if(SolverDataFileName) { |
| // name in .pro file | | // name in .pro file |
|
| if(SolverDataFileName[0] == '/' || SolverDataFileName[0] == '\\'){ | | if(SolverDataFileName[0] == '/' || SolverDataFileName[0] == '\\') { |
| // -> absolute if it starts with '/' or '\' | | // -> absolute if it starts with '/' or '\' |
| strcpy(FileName, SolverDataFileName); | | strcpy(FileName, SolverDataFileName); |
| } | | } |
|
| else{ | | else { |
| // -> relative otherwise FIXME: wrong on Windows - see Fix_RelativePath | | // -> relative otherwise FIXME: wrong on Windows - see Fix_RelativePath |
| strcat(FileName, SolverDataFileName); | | strcat(FileName, SolverDataFileName); |
| } | | } |
| } | | } |
|
| else if (Name_SolverFile){ | | else if(Name_SolverFile) { |
| // name on command line -> always absolute | | // name on command line -> always absolute |
| strcpy(FileName, Name_SolverFile); | | strcpy(FileName, Name_SolverFile); |
| } | | } |
|
| else{ | | else { |
| // default file name -> always relative | | // default file name -> always relative |
| strcat(FileName, Name_DefaultSolverFile); | | strcat(FileName, Name_DefaultSolverFile); |
| } | | } |
| | |
| Message::Info("Loading parameter file '%s'", FileName); | | Message::Info("Loading parameter file '%s'", FileName); |
| | |
|
| init_solver(&Solver->Params, FileName) ; | | init_solver(&Solver->Params, FileName); |
| | |
| i = 0; | | i = 0; |
|
| while(SolverOptions[i] && SolverOptions[i+1]){ | | while(SolverOptions[i] && SolverOptions[i + 1]) { |
| init_solver_option(&Solver->Params, SolverOptions[i], SolverOptions[i+1]) ; | | init_solver_option(&Solver->Params, SolverOptions[i], SolverOptions[i + 1]); |
| i+=2; | | i += 2; |
| } | | } |
| } | | } |
| | |
|
| void LinAlg_SetGlobalSolverOptions(const std::string &opt) | | void LinAlg_SetGlobalSolverOptions(const std::string &opt) {} |
| { | | |
| } | | |
| | |
| void LinAlg_CreateVector(gVector *V, gSolver *Solver, int n) | | void LinAlg_CreateVector(gVector *V, gSolver *Solver, int n) |
| { | | { |
|
| init_vector(n, &V->V) ; | | init_vector(n, &V->V); |
| V->N = n ; | | V->N = n; |
| } | | } |
| | |
| void LinAlg_CreateMatrix(gMatrix *M, gSolver *Solver, int n, int m) | | void LinAlg_CreateMatrix(gMatrix *M, gSolver *Solver, int n, int m) |
| { | | { |
|
| init_matrix(n, &M->M, &Solver->Params) ; | | init_matrix(n, &M->M, &Solver->Params); |
| } | | } |
| | |
| void LinAlg_DestroySolver(gSolver *Solver) | | void LinAlg_DestroySolver(gSolver *Solver) |
| { | | { |
| Message::Debug("'LinAlg_DestroySolver' not yet implemented"); | | Message::Debug("'LinAlg_DestroySolver' not yet implemented"); |
| } | | } |
| | |
|
| void LinAlg_DestroyVector(gVector *V) | | void LinAlg_DestroyVector(gVector *V) { Free(V->V); } |
| { | | |
| Free(V->V); | | |
| } | | |
| | |
|
| void LinAlg_DestroyMatrix(gMatrix *M) | | void LinAlg_DestroyMatrix(gMatrix *M) { free_matrix(&M->M); } |
| { | | |
| free_matrix(&M->M) ; | | |
| } | | |
| | |
|
| void LinAlg_CopyScalar(gScalar *S1, gScalar *S2) | | void LinAlg_CopyScalar(gScalar *S1, gScalar *S2) { S2->s = S1->s; } |
| { | | |
| S2->s = S1->s ; | | |
| } | | |
| | |
| void LinAlg_CopyVector(gVector *V1, gVector *V2) | | void LinAlg_CopyVector(gVector *V1, gVector *V2) |
| { | | { |
|
| memcpy(V2->V, V1->V, V1->N*sizeof(double)) ; | | memcpy(V2->V, V1->V, V1->N * sizeof(double)); |
| } | | } |
| | |
| void LinAlg_SwapVector(gVector *V1, gVector *V2) | | void LinAlg_SwapVector(gVector *V1, gVector *V2) |
| { | | { |
|
| if(V1->N != V2->N){ | | if(V1->N != V2->N) { |
| Message::Error("Cannot swap vectors of different size"); | | Message::Error("Cannot swap vectors of different size"); |
| return; | | return; |
| } | | } |
|
| for(int i = 0; i < V1->N; i++){ | | for(int i = 0; i < V1->N; i++) { |
| double tmp = V1->V[i]; | | double tmp = V1->V[i]; |
| V1->V[i] = V2->V[i]; | | V1->V[i] = V2->V[i]; |
| V2->V[i] = tmp; | | V2->V[i] = tmp; |
| } | | } |
| } | | } |
| | |
| void LinAlg_CopyMatrix(gMatrix *M1, gMatrix *M2) | | void LinAlg_CopyMatrix(gMatrix *M1, gMatrix *M2) |
| { | | { |
| Message::Error("'LinAlg_CopyMatrix' not yet implemented"); | | Message::Error("'LinAlg_CopyMatrix' not yet implemented"); |
| } | | } |
| | |
|
| void LinAlg_ZeroScalar(gScalar *S) | | void LinAlg_ZeroScalar(gScalar *S) { S->s = 0.; } |
| { | | |
| S->s = 0. ; | | |
| } | | |
| | |
|
| void LinAlg_ZeroVector(gVector *V) | | void LinAlg_ZeroVector(gVector *V) { zero_vector(V->N, V->V); } |
| { | | |
| zero_vector(V->N, V->V) ; | | |
| } | | |
| | |
| void LinAlg_ZeroMatrix(gMatrix *M) | | void LinAlg_ZeroMatrix(gMatrix *M) |
| { | | { |
|
| int i ; | | int i; |
| zero_matrix(&M->M) ; | | zero_matrix(&M->M); |
| // la routine de produit matrice vecteur est buggee s'il existe des | | // la routine de produit matrice vecteur est buggee s'il existe des |
| // lignes sans aucun element dans la matrice... | | // 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) ; | | for(i = 0; i < M->M.N; i++) add_matrix_double(&M->M, i + 1, i + 1, 0.0); |
| } | | } |
| | |
|
| void LinAlg_ScanScalar(FILE *file, gScalar *S) | | void LinAlg_ScanScalar(FILE *file, gScalar *S) { fscanf(file, "%lf", &S->s); } |
| { | | |
| fscanf(file, "%lf", &S->s) ; | | |
| } | | |
| | |
| void LinAlg_ScanVector(FILE *file, gVector *V) | | void LinAlg_ScanVector(FILE *file, gVector *V) |
| { | | { |
|
| int i ; | | int i; |
| | |
|
| for(i=0 ; i<V->N ; i++) fscanf(file, "%lf", &V->V[i]) ; | | for(i = 0; i < V->N; i++) fscanf(file, "%lf", &V->V[i]); |
| } | | } |
| | |
| void LinAlg_ScanMatrix(FILE *file, gMatrix *M) | | void LinAlg_ScanMatrix(FILE *file, gMatrix *M) |
| { | | { |
| Message::Error("'LinAlg_ScanMatrix' not yet implemented"); | | Message::Error("'LinAlg_ScanMatrix' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_ReadScalar(FILE *file, gScalar *S) | | void LinAlg_ReadScalar(FILE *file, gScalar *S) |
| { | | { |
| Message::Error("'LinAlg_ReadScalar' not yet implemented"); | | Message::Error("'LinAlg_ReadScalar' not yet implemented"); |
| | |
| skipping to change at line 228 | | skipping to change at line 203 |
|---|
| fread(V->V, sizeof(double), V->N, file); | | fread(V->V, sizeof(double), V->N, file); |
| } | | } |
| | |
| void LinAlg_ReadMatrix(FILE *file, gMatrix *M) | | void LinAlg_ReadMatrix(FILE *file, gMatrix *M) |
| { | | { |
| Message::Error("'LinAlg_ReadMatrix' not yet implemented"); | | Message::Error("'LinAlg_ReadMatrix' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_PrintScalar(FILE *file, gScalar *S) | | void LinAlg_PrintScalar(FILE *file, gScalar *S) |
| { | | { |
|
| fprintf(file, "%.16g", S->s) ; | | fprintf(file, "%.16g", S->s); |
| } | | } |
| | |
| void LinAlg_PrintVector(FILE *file, gVector *V, bool matlab, | | void LinAlg_PrintVector(FILE *file, gVector *V, bool matlab, |
|
| const char* fileName, const char* varName) | | const char *fileName, const char *varName) |
| { | | { |
| if(matlab) Message::Error("Matlab output not available for this vector"); | | if(matlab) Message::Error("Matlab output not available for this vector"); |
|
| formatted_write_vector(file, V->N, V->V, KUL) ; | | formatted_write_vector(file, V->N, V->V, KUL); |
| } | | } |
| | |
| void LinAlg_PrintMatrix(FILE *file, gMatrix *M, bool matlab, | | void LinAlg_PrintMatrix(FILE *file, gMatrix *M, bool matlab, |
|
| const char* fileName, const char* varName) | | const char *fileName, const char *varName) |
| { | | { |
| if(matlab) Message::Error("Matlab output not available for this matrix"); | | if(matlab) Message::Error("Matlab output not available for this matrix"); |
|
| formatted_write_matrix(file, &M->M, KUL) ; | | formatted_write_matrix(file, &M->M, KUL); |
| } | | } |
| | |
| void LinAlg_WriteScalar(FILE *file, gScalar *S) | | void LinAlg_WriteScalar(FILE *file, gScalar *S) |
| { | | { |
| Message::Error("'LinAlg_WriteScalar' not yet implemented"); | | Message::Error("'LinAlg_WriteScalar' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_WriteVector(FILE *file, gVector *V) | | void LinAlg_WriteVector(FILE *file, gVector *V) |
| { | | { |
| fwrite(V->V, sizeof(double), V->N, file); | | fwrite(V->V, sizeof(double), V->N, file); |
| fprintf(file, "\n"); | | fprintf(file, "\n"); |
| } | | } |
| | |
| void LinAlg_WriteMatrix(FILE *file, gMatrix *M) | | void LinAlg_WriteMatrix(FILE *file, gMatrix *M) |
| { | | { |
|
| binary_write_matrix(&M->M, "A", ".mat") ; | | binary_write_matrix(&M->M, "A", ".mat"); |
| } | | } |
| | |
|
| void LinAlg_GetVectorSize(gVector *V, int *i) | | void LinAlg_GetVectorSize(gVector *V, int *i) { *i = V->N; } |
| { | | |
| *i = V->N ; | | |
| } | | |
| | |
| void LinAlg_GetLocalVectorRange(gVector *V, int *low, int *high) | | void LinAlg_GetLocalVectorRange(gVector *V, int *low, int *high) |
| { | | { |
|
| *low = 0 ; | | *low = 0; |
| *high = V->N ; | | *high = V->N; |
| } | | } |
| | |
|
| void LinAlg_GetMatrixSize(gMatrix *M, int *i, int *j) | | void LinAlg_GetMatrixSize(gMatrix *M, int *i, int *j) { *i = *j = M->M.N; } |
| { | | |
| *i = *j = M->M.N ; | | |
| } | | |
| | |
| void LinAlg_GetLocalMatrixRange(gMatrix *M, int *low, int *high) | | void LinAlg_GetLocalMatrixRange(gMatrix *M, int *low, int *high) |
| { | | { |
|
| *low = 0 ; | | *low = 0; |
| *high = M->M.N ; | | *high = M->M.N; |
| } | | } |
| | |
|
| void LinAlg_GetDoubleInScalar(double *d, gScalar *S) | | void LinAlg_GetDoubleInScalar(double *d, gScalar *S) { *d = S->s; } |
| { | | |
| *d = S->s ; | | |
| } | | |
| | |
| void LinAlg_GetComplexInScalar(double *d1, double *d2, gScalar *S) | | void LinAlg_GetComplexInScalar(double *d1, double *d2, gScalar *S) |
| { | | { |
| Message::Error("'LinAlg_GetComplexInScalar' not available with this Solver"); | | Message::Error("'LinAlg_GetComplexInScalar' not available with this Solver"); |
| } | | } |
| | |
|
| void LinAlg_GetScalarInVector(gScalar *S, gVector *V, int i) | | void LinAlg_GetScalarInVector(gScalar *S, gVector *V, int i) { S->s = V->V[i]; } |
| { | | |
| S->s = V->V[i] ; | | |
| } | | |
| | |
|
| void LinAlg_GetDoubleInVector(double *d, gVector *V, int i) | | void LinAlg_GetDoubleInVector(double *d, gVector *V, int i) { *d = V->V[i]; } |
| { | | |
| *d = V->V[i] ; | | |
| } | | |
| | |
| void LinAlg_GetAbsDoubleInVector(double *d, gVector *V, int i) | | void LinAlg_GetAbsDoubleInVector(double *d, gVector *V, int i) |
| { | | { |
|
| *d = fabs(V->V[i]) ; | | *d = fabs(V->V[i]); |
| } | | } |
| | |
| void LinAlg_GetComplexInVector(double *d1, double *d2, gVector *V, int i, int j) | | void LinAlg_GetComplexInVector(double *d1, double *d2, gVector *V, int i, int j) |
| { | | { |
|
| *d1 = V->V[i] ; | | *d1 = V->V[i]; |
| *d2 = V->V[j] ; | | *d2 = V->V[j]; |
| } | | } |
| | |
| void LinAlg_GetScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) | | void LinAlg_GetScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) |
| { | | { |
| Message::Error("'LinAlg_GetScalarInMatrix' not yet implemented"); | | Message::Error("'LinAlg_GetScalarInMatrix' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_GetDoubleInMatrix(double *d, gMatrix *M, int i, int j) | | void LinAlg_GetDoubleInMatrix(double *d, gMatrix *M, int i, int j) |
| { | | { |
| get_element_in_matrix(&M->M, i, j, d); | | get_element_in_matrix(&M->M, i, j, d); |
| } | | } |
| | |
|
| void LinAlg_GetComplexInMatrix(double *d1, double *d2, gMatrix *M, int i, int j, | | void LinAlg_GetComplexInMatrix(double *d1, double *d2, gMatrix *M, int i, int j, |
| int k, int l) | | int k, int l) |
| { | | { |
| Message::Error("'LinAlg_GetComplexInMatrix' not yet implemented"); | | Message::Error("'LinAlg_GetComplexInMatrix' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_GetColumnInMatrix(gMatrix *M, int col, gVector *V1) | | void LinAlg_GetColumnInMatrix(gMatrix *M, int col, gVector *V1) |
| { | | { |
| get_column_in_matrix(&M->M, col, V1->V); | | get_column_in_matrix(&M->M, col, V1->V); |
| } | | } |
| | |
|
| void LinAlg_SetScalar(gScalar *S, double *d) | | void LinAlg_SetScalar(gScalar *S, double *d) { S->s = d[0]; } |
| { | | |
| S->s = d[0] ; | | |
| } | | |
| | |
| void LinAlg_SetVector(gVector *V, double *v) | | void LinAlg_SetVector(gVector *V, double *v) |
| { | | { |
| int i; | | int i; |
| | |
|
| for(i=0; i<V->N; i++) V->V[i] = *v ; | | for(i = 0; i < V->N; i++) V->V[i] = *v; |
| } | | } |
| | |
|
| void LinAlg_SetScalarInVector(gScalar *S, gVector *V, int i) | | void LinAlg_SetScalarInVector(gScalar *S, gVector *V, int i) { V->V[i] = S->s; } |
| { | | |
| V->V[i] = S->s ; | | |
| } | | |
| | |
|
| void LinAlg_SetDoubleInVector(double d, gVector *V, int i) | | void LinAlg_SetDoubleInVector(double d, gVector *V, int i) { V->V[i] = d; } |
| { | | |
| V->V[i] = d ; | | |
| } | | |
| | |
| void LinAlg_SetComplexInVector(double d1, double d2, gVector *V, int i, int j) | | void LinAlg_SetComplexInVector(double d1, double d2, gVector *V, int i, int j) |
| { | | { |
|
| V->V[i] = d1 ; | | V->V[i] = d1; |
| V->V[j] = d2 ; | | V->V[j] = d2; |
| } | | } |
| | |
| void LinAlg_SetScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) | | void LinAlg_SetScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) |
| { | | { |
| Message::Error("'LinAlg_SetScalarInMatrix' not yet implemented"); | | Message::Error("'LinAlg_SetScalarInMatrix' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_SetDoubleInMatrix(double d, gMatrix *M, int i, int j) | | void LinAlg_SetDoubleInMatrix(double d, gMatrix *M, int i, int j) |
| { | | { |
| Message::Error("'LinAlg_SetDoubleInMatrix' not yet implemented"); | | Message::Error("'LinAlg_SetDoubleInMatrix' not yet implemented"); |
| } | | } |
| | |
|
| void LinAlg_SetComplexInMatrix(double d1, double d2, gMatrix *M, int i, int j, i | | void LinAlg_SetComplexInMatrix(double d1, double d2, gMatrix *M, int i, int j, |
| nt k, int l) | | int k, int l) |
| { | | { |
| Message::Error("'LinAlg_SetComplexInMatrix' not yet implemented"); | | Message::Error("'LinAlg_SetComplexInMatrix' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_AddScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) | | void LinAlg_AddScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) |
| { | | { |
|
| S3->s = S1->s + S2->s ; | | S3->s = S1->s + S2->s; |
| } | | } |
| | |
| void LinAlg_DummyVector(gVector *V) | | void LinAlg_DummyVector(gVector *V) |
| { | | { |
|
| int * DummyDof, i; | | int *DummyDof, i; |
| | |
| DummyDof = Current.DofData->DummyDof; | | DummyDof = Current.DofData->DummyDof; |
|
| if (DummyDof == NULL) return ; | | if(DummyDof == NULL) return; |
| | |
|
| for (i=0 ; i<V->N ; i++) if (DummyDof[i] == 1) V->V[i] = 0 ; | | for(i = 0; i < V->N; i++) |
| | if(DummyDof[i] == 1) V->V[i] = 0; |
| } | | } |
| | |
| void LinAlg_AddScalarInVector(gScalar *S, gVector *V, int i) | | void LinAlg_AddScalarInVector(gScalar *S, gVector *V, int i) |
| { | | { |
|
| int * DummyDof; | | int *DummyDof; |
| | |
|
| if ((DummyDof = Current.DofData->DummyDof)) | | if((DummyDof = Current.DofData->DummyDof)) |
| if (DummyDof[i] == 1) return ; | | if(DummyDof[i] == 1) return; |
| | |
|
| V->V[i] += S->s ; | | V->V[i] += S->s; |
| } | | } |
| | |
| void LinAlg_AddDoubleInVector(double d, gVector *V, int i) | | void LinAlg_AddDoubleInVector(double d, gVector *V, int i) |
| { | | { |
|
| int * DummyDof; | | int *DummyDof; |
| | |
|
| if ((DummyDof = Current.DofData->DummyDof)) | | if((DummyDof = Current.DofData->DummyDof)) |
| if (DummyDof[i] == 1) | | if(DummyDof[i] == 1) return; |
| return ; | | |
| | |
|
| V->V[i] += d ; | | V->V[i] += d; |
| } | | } |
| | |
| void LinAlg_AddComplexInVector(double d1, double d2, gVector *V, int i, int j) | | void LinAlg_AddComplexInVector(double d1, double d2, gVector *V, int i, int j) |
| { | | { |
|
| int * DummyDof, iok,jok; | | int *DummyDof, iok, jok; |
| | |
|
| iok=jok=1; | | iok = jok = 1; |
| | |
|
| if ((DummyDof = Current.DofData->DummyDof)) { | | if((DummyDof = Current.DofData->DummyDof)) { |
| if (DummyDof[i] == 1) iok=0; | | if(DummyDof[i] == 1) iok = 0; |
| if (DummyDof[j] == 1) jok=0; | | if(DummyDof[j] == 1) jok = 0; |
| } | | } |
| | |
|
| if (iok) | | if(iok) V->V[i] += d1; |
| V->V[i] += d1 ; | | |
| | |
|
| if (jok) | | if(jok) V->V[j] += d2; |
| V->V[j] += d2 ; | | |
| } | | } |
| | |
| void LinAlg_AddScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) | | void LinAlg_AddScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) |
| { | | { |
|
| | int *DummyDof; |
| | |
|
| int * DummyDof; | | if((DummyDof = Current.DofData->DummyDof)) |
| | if((DummyDof[i] == 1 || DummyDof[j] == 1) && (i != j)) return; |
| | |
|
| if ((DummyDof = Current.DofData->DummyDof)) | | add_matrix_double(&M->M, i + 1, j + 1, S->s); |
| if ( (DummyDof[i] == 1 || DummyDof[j] == 1) && (i != j) ) | | |
| return ; | | |
| | | |
| add_matrix_double(&M->M, i+1, j+1, S->s) ; | | |
| } | | } |
| | |
| void LinAlg_AddDoubleInMatrix(double d, gMatrix *M, int i, int j) | | void LinAlg_AddDoubleInMatrix(double d, gMatrix *M, int i, int j) |
| { | | { |
|
| int * DummyDof; | | int *DummyDof; |
| | |
|
| if ((DummyDof = Current.DofData->DummyDof)) | | if((DummyDof = Current.DofData->DummyDof)) |
| if ( (DummyDof[i] == 1 || DummyDof[j] == 1) && (i != j) ) | | if((DummyDof[i] == 1 || DummyDof[j] == 1) && (i != j)) return; |
| return ; | | |
| | |
|
| add_matrix_double(&M->M, i+1, j+1, d) ; | | add_matrix_double(&M->M, i + 1, j + 1, d); |
| } | | } |
| | |
|
| void LinAlg_AddComplexInMatrix(double d1, double d2, gMatrix *M, int i, int j, i | | void LinAlg_AddComplexInMatrix(double d1, double d2, gMatrix *M, int i, int j, |
| nt k, int l) | | int k, int l) |
| { | | { |
|
| if(d1){ | | if(d1) { |
| add_matrix_double(&M->M, i+1, j+1, d1) ; | | add_matrix_double(&M->M, i + 1, j + 1, d1); |
| add_matrix_double(&M->M, k+1, l+1, d1) ; | | add_matrix_double(&M->M, k + 1, l + 1, d1); |
| } | | } |
|
| if(d2){ | | if(d2) { |
| add_matrix_double(&M->M, i+1, l+1, -d2) ; | | add_matrix_double(&M->M, i + 1, l + 1, -d2); |
| add_matrix_double(&M->M, k+1, j+1, d2) ; | | add_matrix_double(&M->M, k + 1, j + 1, d2); |
| } | | } |
| } | | } |
| | |
| void LinAlg_AddVectorVector(gVector *V1, gVector *V2, gVector *V3) | | void LinAlg_AddVectorVector(gVector *V1, gVector *V2, gVector *V3) |
| { | | { |
| if(V3 == V1) | | if(V3 == V1) |
|
| add_vector_vector(V1->N, V1->V, V2->V) ; | | add_vector_vector(V1->N, V1->V, V2->V); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_AddVectorVector'"); | | Message::Error("Wrong arguments in 'LinAlg_AddVectorVector'"); |
| } | | } |
| | |
|
| void LinAlg_AddVectorProdVectorDouble(gVector *V1, gVector *V2, double d, gVecto | | void LinAlg_AddVectorProdVectorDouble(gVector *V1, gVector *V2, double d, |
| r *V3) | | gVector *V3) |
| { | | { |
| if(V3 == V1) | | if(V3 == V1) |
|
| add_vector_prod_vector_double(V1->N, V1->V, V2->V, d) ; | | add_vector_prod_vector_double(V1->N, V1->V, V2->V, d); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_AddVectorProdVectorDouble'"); | | Message::Error("Wrong arguments in 'LinAlg_AddVectorProdVectorDouble'"); |
| } | | } |
| | |
| void LinAlg_AddProdVectorDoubleProdVectorDouble(double alpha, gVector *V1, | | void LinAlg_AddProdVectorDoubleProdVectorDouble(double alpha, gVector *V1, |
|
| double beta, gVector *V2, gVector *V3) | | double beta, gVector *V2, |
| | gVector *V3) |
| { | | { |
|
| Message::Error("'LinAlg_AddProdVectorDoubleProdVectorDouble' not yet implement | | Message::Error( |
| ed"); | | "'LinAlg_AddProdVectorDoubleProdVectorDouble' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_AddMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3) | | void LinAlg_AddMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3) |
| { | | { |
| if(M3 == M1) | | if(M3 == M1) |
|
| add_matrix_matrix(&M1->M, &M2->M) ; | | add_matrix_matrix(&M1->M, &M2->M); |
| else if(M3 == M2) | | else if(M3 == M2) |
|
| add_matrix_matrix(&M2->M, &M1->M) ; | | add_matrix_matrix(&M2->M, &M1->M); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_AddMatrixMatrix'"); | | Message::Error("Wrong arguments in 'LinAlg_AddMatrixMatrix'"); |
| } | | } |
| | |
|
| void LinAlg_AddMatrixProdMatrixDouble(gMatrix *M1, gMatrix *M2, double d, gMatri | | void LinAlg_AddMatrixProdMatrixDouble(gMatrix *M1, gMatrix *M2, double d, |
| x *M3) | | gMatrix *M3) |
| { | | { |
| if(M3 == M1) | | if(M3 == M1) |
|
| add_matrix_prod_matrix_double(&M1->M, &M2->M, d) ; | | add_matrix_prod_matrix_double(&M1->M, &M2->M, d); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_AddMatrixProdMatrixDouble'"); | | Message::Error("Wrong arguments in 'LinAlg_AddMatrixProdMatrixDouble'"); |
| } | | } |
| | |
| void LinAlg_SubScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) | | void LinAlg_SubScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) |
| { | | { |
|
| S3->s = S1->s - S2->s ; | | S3->s = S1->s - S2->s; |
| } | | } |
| | |
| void LinAlg_SubVectorVector(gVector *V1, gVector *V2, gVector *V3) | | void LinAlg_SubVectorVector(gVector *V1, gVector *V2, gVector *V3) |
| { | | { |
| if(V3 == V1) | | if(V3 == V1) |
|
| sub_vector_vector_1(V1->N, V1->V, V2->V) ; | | sub_vector_vector_1(V1->N, V1->V, V2->V); |
| else if (V3 == V2) | | else if(V3 == V2) |
| sub_vector_vector_2(V1->N, V1->V, V2->V) ; | | sub_vector_vector_2(V1->N, V1->V, V2->V); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_SubVectorVector'"); | | Message::Error("Wrong arguments in 'LinAlg_SubVectorVector'"); |
| } | | } |
| | |
| void LinAlg_SubMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3) | | void LinAlg_SubMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3) |
| { | | { |
| Message::Error("'LinAlg_SubMatrixMatrix' not yet implemented"); | | Message::Error("'LinAlg_SubMatrixMatrix' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_ProdScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) | | void LinAlg_ProdScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) |
| { | | { |
|
| S3->s = S1->s * S2->s ; | | S3->s = S1->s * S2->s; |
| } | | } |
| | |
| void LinAlg_ProdScalarDouble(gScalar *S1, double d, gScalar *S2) | | void LinAlg_ProdScalarDouble(gScalar *S1, double d, gScalar *S2) |
| { | | { |
|
| S2->s = S1->s * d ; | | S2->s = S1->s * d; |
| } | | } |
| | |
|
| void LinAlg_ProdScalarComplex(gScalar *S, double d1, double d2, double *d3, doub | | void LinAlg_ProdScalarComplex(gScalar *S, double d1, double d2, double *d3, |
| le *d4) | | double *d4) |
| { | | { |
|
| *d3 = S->s * d1 ; | | *d3 = S->s * d1; |
| *d4 = S->s * d2 ; | | *d4 = S->s * d2; |
| } | | } |
| | |
| void LinAlg_ProdVectorScalar(gVector *V1, gScalar *S, gVector *V2) | | void LinAlg_ProdVectorScalar(gVector *V1, gScalar *S, gVector *V2) |
| { | | { |
| Message::Error("'LinAlg_ProdVectorScalar' not yet implemented"); | | Message::Error("'LinAlg_ProdVectorScalar' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_ProdVectorDouble(gVector *V1, double d, gVector *V2) | | void LinAlg_ProdVectorDouble(gVector *V1, double d, gVector *V2) |
| { | | { |
| if(V2 == V1) | | if(V2 == V1) |
| | |
| skipping to change at line 562 | | skipping to change at line 516 |
|---|
| Message::Error("Wrong arguments in 'LinAlg_ProdVectorDouble'"); | | Message::Error("Wrong arguments in 'LinAlg_ProdVectorDouble'"); |
| } | | } |
| | |
| void LinAlg_ProdVectorComplex(gVector *V1, double d1, double d2, gVector *V2) | | void LinAlg_ProdVectorComplex(gVector *V1, double d1, double d2, gVector *V2) |
| { | | { |
| Message::Error("'LinAlg_ProdVectorComplex' not yet implemented"); | | Message::Error("'LinAlg_ProdVectorComplex' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_ProdVectorVector(gVector *V1, gVector *V2, double *d) | | void LinAlg_ProdVectorVector(gVector *V1, gVector *V2, double *d) |
| { | | { |
|
| prodsc_vector_vector (V1->N, V1->V, V2->V, d) ; | | prodsc_vector_vector(V1->N, V1->V, V2->V, d); |
| } | | } |
| | |
| void LinAlg_ProdMatrixVector(gMatrix *M, gVector *V1, gVector *V2) | | void LinAlg_ProdMatrixVector(gMatrix *M, gVector *V1, gVector *V2) |
| { | | { |
| if(V2 == V1) | | if(V2 == V1) |
| Message::Error("Wrong arguments in 'LinAlg_ProdMatrixVector'"); | | Message::Error("Wrong arguments in 'LinAlg_ProdMatrixVector'"); |
| else | | else |
| prod_matrix_vector(&M->M, V1->V, V2->V); | | prod_matrix_vector(&M->M, V1->V, V2->V); |
| } | | } |
| | |
| void LinAlg_ProdMatrixScalar(gMatrix *M1, gScalar *S, gMatrix *M2) | | void LinAlg_ProdMatrixScalar(gMatrix *M1, gScalar *S, gMatrix *M2) |
| { | | { |
| if(M2 == M1) | | if(M2 == M1) |
|
| prod_matrix_double (&M1->M, S->s); | | prod_matrix_double(&M1->M, S->s); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_ProdMatrixScalar'"); | | Message::Error("Wrong arguments in 'LinAlg_ProdMatrixScalar'"); |
| } | | } |
| | |
| void LinAlg_ProdMatrixDouble(gMatrix *M1, double d, gMatrix *M2) | | void LinAlg_ProdMatrixDouble(gMatrix *M1, double d, gMatrix *M2) |
| { | | { |
| if(M2 == M1) | | if(M2 == M1) |
|
| prod_matrix_double (&M1->M, d); | | prod_matrix_double(&M1->M, d); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_ProdMatrixDouble'"); | | Message::Error("Wrong arguments in 'LinAlg_ProdMatrixDouble'"); |
| } | | } |
| | |
| void LinAlg_ProdMatrixComplex(gMatrix *M1, double d1, double d2, gMatrix *M2) | | void LinAlg_ProdMatrixComplex(gMatrix *M1, double d1, double d2, gMatrix *M2) |
| { | | { |
| Message::Error("'LinAlg_ProdMatrixComplex' not yet implemented"); | | Message::Error("'LinAlg_ProdMatrixComplex' not yet implemented"); |
| } | | } |
| | |
| void LinAlg_DivScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) | | void LinAlg_DivScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) |
| { | | { |
|
| S3->s = S1->s / S2->s ; | | S3->s = S1->s / S2->s; |
| } | | } |
| | |
| void LinAlg_DivScalarDouble(gScalar *S1, double d, gScalar *S2) | | void LinAlg_DivScalarDouble(gScalar *S1, double d, gScalar *S2) |
| { | | { |
|
| S2->s = S1->s / d ; | | S2->s = S1->s / d; |
| } | | } |
| | |
| void LinAlg_VectorNorm2(gVector *V1, double *norm) | | void LinAlg_VectorNorm2(gVector *V1, double *norm) |
| { | | { |
| norm2_vector(V1->N, V1->V, norm); | | norm2_vector(V1->N, V1->V, norm); |
| } | | } |
| | |
| void LinAlg_VectorNormInf(gVector *V1, double *norm) | | void LinAlg_VectorNormInf(gVector *V1, double *norm) |
| { | | { |
| norminf_vector(V1->N, V1->V, norm); | | norminf_vector(V1->N, V1->V, norm); |
| } | | } |
| | |
|
| void LinAlg_AssembleMatrix(gMatrix *M) | | void LinAlg_AssembleMatrix(gMatrix *M) {} |
| { | | |
| } | | |
| | |
|
| void LinAlg_AssembleVector(gVector *V) | | void LinAlg_AssembleVector(gVector *V) {} |
| { | | |
| } | | |
| | |
|
| void LinAlg_Solve(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, int solve | | void LinAlg_Solve(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, |
| rIndex) | | int solverIndex) |
| { | | { |
| solve_matrix(&A->M, &Solver->Params, B->V, X->V); | | solve_matrix(&A->M, &Solver->Params, B->V, X->V); |
| } | | } |
| | |
|
| void LinAlg_SolveAgain(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, int | | void LinAlg_SolveAgain(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, |
| solverIndex) | | int solverIndex) |
| { | | { |
| int tmp = Solver->Params.Re_Use_LU; | | int tmp = Solver->Params.Re_Use_LU; |
| Solver->Params.Re_Use_LU = 1; | | Solver->Params.Re_Use_LU = 1; |
| solve_matrix(&A->M, &Solver->Params, B->V, X->V); | | solve_matrix(&A->M, &Solver->Params, B->V, X->V); |
| Solver->Params.Re_Use_LU = tmp; | | Solver->Params.Re_Use_LU = tmp; |
| } | | } |
| | |
|
| void LinAlg_SolveNL(gMatrix *A, gVector *B, gMatrix *J, gVector *R, gSolver *Sol | | void LinAlg_SolveNL(gMatrix *A, gVector *B, gMatrix *J, gVector *R, |
| ver, | | gSolver *Solver, gVector *X, int solverIndex) |
| gVector *X, int solverIndex) | | |
| { | | { |
| Message::Error("'LinAlg_SolveNL' not yet implemented for Sparskit"); | | Message::Error("'LinAlg_SolveNL' not yet implemented for Sparskit"); |
| } | | } |
| | |
| #endif | | #endif |
| | | |
| End of changes. 108 change blocks. |
|---|
| 219 lines changed or deleted | | 160 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 