"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Kernel/EigenSolve_ARPACK.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

[ To the main getdp source changes report ]

EigenSolve_ARPACK.cpp  (getdp-3.4.0-source.tgz):EigenSolve_ARPACK.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):
// Alexandru Mustatea // Alexandru Mustatea
// Andre Nicolet // Andre Nicolet
// //
#include "GetDPConfig.h" #include "GetDPConfig.h"
skipping to change at line 24 skipping to change at line 24
#include <string.h> #include <string.h>
#include <math.h> #include <math.h>
#include "Message.h" #include "Message.h"
#include "ProData.h" #include "ProData.h"
#include "DofData.h" #include "DofData.h"
#include "Cal_Quantity.h" #include "Cal_Quantity.h"
#include "MallocUtils.h" #include "MallocUtils.h"
#include "OS.h" #include "OS.h"
#define SQU(a) ((a)*(a)) #define SQU(a) ((a) * (a))
#define TWO_PI 6.2831853071795865 #define TWO_PI 6.2831853071795865
extern struct CurrentData Current ; extern struct CurrentData Current;
extern char *Name_Path ; extern char *Name_Path;
struct EigenPar { struct EigenPar {
double prec; double prec;
int size; int size;
int reortho; int reortho;
} ; };
#if defined(HAVE_NO_UNDERSCORE) #if defined(HAVE_NO_UNDERSCORE)
#define znaupd_ znaupd #define znaupd_ znaupd
#define zneupd_ zneupd #define zneupd_ zneupd
#endif #endif
extern "C" { extern "C" {
typedef struct {double re; double im;} complex_16; typedef struct {
extern void znaupd_(int *ido, char *bmat, int *n, double re;
char *which, int *nev, double *tol, complex_16 resid[], double im;
int *ncv, complex_16 v[], int *ldv, int iparam[], } complex_16;
int ipntr[], complex_16 workd[], complex_16 workl[], int *l extern void znaupd_(int *ido, char *bmat, int *n, char *which, int *nev,
workl, double *tol, complex_16 resid[], int *ncv, complex_16 v[],
double rwork[], int *info); int *ldv, int iparam[], int ipntr[], complex_16 workd[],
extern void zneupd_(unsigned int *rvec, char *howmny, unsigned int select[], complex_16 workl[], int *lworkl, double rwork[], int *info);
complex_16 d[], complex_16 z[], int *ldz, complex_16 *sigma extern void zneupd_(unsigned int *rvec, char *howmny, unsigned int select[],
, complex_16 d[], complex_16 z[], int *ldz, complex_16 *sigma,
complex_16 workev[], char *bmat, int *n, char *which, complex_16 workev[], char *bmat, int *n, char *which,
int *nev, double *tol, complex_16 resid[], int *ncv, int *nev, double *tol, complex_16 resid[], int *ncv,
complex_16 v[], int *ldv, int iparam[], int ipntr[], complex_16 v[], int *ldv, int iparam[], int ipntr[],
complex_16 workd[], complex_16 workl[], int *lworkl, complex_16 workd[], complex_16 workl[], int *lworkl,
double rwork[], int *info); double rwork[], int *info);
} }
static void EigenGetDouble(const char *text, double *d) static void EigenGetDouble(const char *text, double *d)
{ {
char str[256]; char str[256];
printf("%s (default=%.16g): ", text, *d); printf("%s (default=%.16g): ", text, *d);
if(fgets(str, sizeof(str), stdin)){ if(fgets(str, sizeof(str), stdin)) {
if(strlen(str) && strcmp(str, "\n")) if(strlen(str) && strcmp(str, "\n")) *d = atof(str);
*d = atof(str);
} }
} }
static void EigenGetInt(const char *text, int *i) static void EigenGetInt(const char *text, int *i)
{ {
char str[256]; char str[256];
printf("%s (default=%d): ", text, *i); printf("%s (default=%d): ", text, *i);
if(fgets(str, sizeof(str), stdin)){ if(fgets(str, sizeof(str), stdin)) {
if(strlen(str) && strcmp(str, "\n")) if(strlen(str) && strcmp(str, "\n")) *i = atoi(str);
*i = atoi(str);
} }
} }
void EigenPar(const char *filename, struct EigenPar *par) void EigenPar(const char *filename, struct EigenPar *par)
{ {
char path[1024]; char path[1024];
FILE *fp; FILE *fp;
/* set some defaults */ /* set some defaults */
par->prec = 1.e-4; par->prec = 1.e-4;
par->reortho = 0; par->reortho = 0;
par->size = 50; par->size = 50;
/* try to read parameters from file */ /* try to read parameters from file */
strcpy(path, Name_Path); strcpy(path, Name_Path);
strcat(path, filename); strcat(path, filename);
fp = FOpen(path, "r"); fp = FOpen(path, "r");
if(fp) { if(fp) {
Message::Info("Loading eigenproblem parameter file '%s'", path); Message::Info("Loading eigenproblem parameter file '%s'", path);
if(fscanf(fp, "%lf %d %d", &par->prec, &par->reortho, &par->size) != 3){ if(fscanf(fp, "%lf %d %d", &par->prec, &par->reortho, &par->size) != 3) {
Message::Error("Could not read parameters"); Message::Error("Could not read parameters");
} }
fclose(fp); fclose(fp);
} }
else{ else {
fp = FOpen(path, "w"); fp = FOpen(path, "w");
if(fp){ if(fp) {
if(!Message::UseOnelab()){ if(!Message::UseOnelab()) {
/* get parameters from command line */ /* get parameters from command line */
EigenGetDouble("Precision", &par->prec); EigenGetDouble("Precision", &par->prec);
EigenGetInt("Reorthogonalization", &par->reortho); EigenGetInt("Reorthogonalization", &par->reortho);
EigenGetInt("Krylov basis size", &par->size); EigenGetInt("Krylov basis size", &par->size);
} }
/* write file */ /* write file */
fprintf(fp, "%.16g\n", par->prec); fprintf(fp, "%.16g\n", par->prec);
fprintf(fp, "%d\n", par->reortho); fprintf(fp, "%d\n", par->reortho);
fprintf(fp, "%d\n", par->size); fprintf(fp, "%d\n", par->size);
fprintf(fp, fprintf(
"/*\n" fp, "/*\n"
" The numbers above are the parameters for the numerical\n" " The numbers above are the parameters for the numerical\n"
" eigenvalue problem:\n" " eigenvalue problem:\n"
"\n" "\n"
" prec = aimed accuracy for eigenvectors (default=1.e-4)\n" " prec = aimed accuracy for eigenvectors (default=1.e-4)\n"
" reortho = reorthogonalisation of Krylov basis: yes=1, no=0 (def " reortho = reorthogonalisation of Krylov basis: yes=1, no=0 "
ault=0) \n" "(default=0) \n"
" size = size of the Krylov basis\n" " size = size of the Krylov basis\n"
"\n" "\n"
" The shift is given in the .pro file because its choice relies\n " The shift is given in the .pro file because its choice relies\n"
" " on physical considerations.\n"
" on physical considerations.\n" "*/");
"*/");
fclose(fp); fclose(fp);
} }
else{ else {
Message::Error("Unable to open file '%s'", path); Message::Error("Unable to open file '%s'", path);
} }
} }
Message::Info("Eigenproblem parameters: prec = %g, reortho = %d, size = %d", Message::Info("Eigenproblem parameters: prec = %g, reortho = %d, size = %d",
par->prec, par->reortho, par->size); par->prec, par->reortho, par->size);
} }
/* This routine uses Arpack to solve Generalized Complex Non-Hermitian /* This routine uses Arpack to solve Generalized Complex Non-Hermitian
eigenvalue problems. We don't use the "Generalized" Arpack mode eigenvalue problems. We don't use the "Generalized" Arpack mode
(bmat=='G') since it requires M to be Hermitian. Instead, we use (bmat=='G') since it requires M to be Hermitian. Instead, we use
the regular mode (bmat='I') and apply the shift "by hand", which the regular mode (bmat='I') and apply the shift "by hand", which
allows us to use arbitrary matrices K and M. */ allows us to use arbitrary matrices K and M. */
static void Arpack2GetDP(int N, complex_16 *in, gVector *out) static void Arpack2GetDP(int N, complex_16 *in, gVector *out)
{ {
int i, j; int i, j;
double re, im; double re, im;
int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1; int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1;
for(i = 0; i < N; i++){ for(i = 0; i < N; i++) {
re = in[i].re; re = in[i].re;
im = in[i].im; im = in[i].im;
j = i * incr; j = i * incr;
if(Current.NbrHar == 2) if(Current.NbrHar == 2)
LinAlg_SetComplexInVector(re, im, out, j, j+1); LinAlg_SetComplexInVector(re, im, out, j, j + 1);
else else
LinAlg_SetDoubleInVector(re, out, j); LinAlg_SetDoubleInVector(re, out, j);
} }
LinAlg_AssembleVector(out); LinAlg_AssembleVector(out);
} }
static void Arpack2GetDPSplit(int N, complex_16 *in, gVector *out1, gVector *out static void Arpack2GetDPSplit(int N, complex_16 *in, gVector *out1,
2) gVector *out2)
{ {
int i, j; int i, j;
double re, im; double re, im;
int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1; int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1;
for(i = 0; i < N/2; i++){ for(i = 0; i < N / 2; i++) {
j = i * incr; j = i * incr;
re = in[i].re; re = in[i].re;
im = in[i].im; im = in[i].im;
if(Current.NbrHar == 2) if(Current.NbrHar == 2)
LinAlg_SetComplexInVector(re, im, out1, j, j+1); LinAlg_SetComplexInVector(re, im, out1, j, j + 1);
else else
LinAlg_SetDoubleInVector(re, out1, j); LinAlg_SetDoubleInVector(re, out1, j);
re = in[N/2+i].re; re = in[N / 2 + i].re;
im = in[N/2+i].im; im = in[N / 2 + i].im;
if(Current.NbrHar == 2) if(Current.NbrHar == 2)
LinAlg_SetComplexInVector(re, im, out2, j, j+1); LinAlg_SetComplexInVector(re, im, out2, j, j + 1);
else else
LinAlg_SetDoubleInVector(re, out2, j); LinAlg_SetDoubleInVector(re, out2, j);
} }
LinAlg_AssembleVector(out1); LinAlg_AssembleVector(out1);
LinAlg_AssembleVector(out2); LinAlg_AssembleVector(out2);
} }
static void GetDP2Arpack(gVector *in, complex_16 *out) static void GetDP2Arpack(gVector *in, complex_16 *out)
{ {
int i, N; int i, N;
double re, im = 0.; double re, im = 0.;
int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1; int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1;
LinAlg_GetVectorSize(in, &N); LinAlg_GetVectorSize(in, &N);
for(i = 0; i < N; i += incr){ for(i = 0; i < N; i += incr) {
if(Current.NbrHar == 2) if(Current.NbrHar == 2)
LinAlg_GetComplexInVector(&re, &im, in, i, i+1); LinAlg_GetComplexInVector(&re, &im, in, i, i + 1);
else else
LinAlg_GetDoubleInVector(&re, in, i); LinAlg_GetDoubleInVector(&re, in, i);
out[i/incr].re = re; out[i / incr].re = re;
out[i/incr].im = im; out[i / incr].im = im;
} }
} }
static void GetDP2ArpackMerge(gVector *in1, gVector *in2, complex_16 *out) static void GetDP2ArpackMerge(gVector *in1, gVector *in2, complex_16 *out)
{ {
int i, N; int i, N;
double re, im = 0.; double re, im = 0.;
int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1; int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1;
LinAlg_GetVectorSize(in1, &N); LinAlg_GetVectorSize(in1, &N);
for(i = 0; i < N; i += incr){ for(i = 0; i < N; i += incr) {
if(Current.NbrHar == 2) if(Current.NbrHar == 2)
LinAlg_GetComplexInVector(&re, &im, in1, i, i+1); LinAlg_GetComplexInVector(&re, &im, in1, i, i + 1);
else else
LinAlg_GetDoubleInVector(&re, in1, i); LinAlg_GetDoubleInVector(&re, in1, i);
out[i/incr].re = re; out[i / incr].re = re;
out[i/incr].im = im; out[i / incr].im = im;
if(Current.NbrHar == 2) if(Current.NbrHar == 2)
LinAlg_GetComplexInVector(&re, &im, in2, i, i+1); LinAlg_GetComplexInVector(&re, &im, in2, i, i + 1);
else else
LinAlg_GetDoubleInVector(&re, in2, i); LinAlg_GetDoubleInVector(&re, in2, i);
out[N/incr + i/incr].re = re; out[N / incr + i / incr].re = re;
out[N/incr + i/incr].im = im; out[N / incr + i / incr].im = im;
} }
} }
void EigenSolve_ARPACK(struct DofData * DofData_P, int NumEigenvalues, void EigenSolve_ARPACK(struct DofData *DofData_P, int NumEigenvalues,
double shift_r, double shift_i, int FilterExpressionIndex double shift_r, double shift_i,
) int FilterExpressionIndex)
{ {
struct EigenPar eigenpar; struct EigenPar eigenpar;
struct Solution Solution_S; struct Solution Solution_S;
gVector v1, v2, w1, w2, x, y; gVector v1, v2, w1, w2, x, y;
int n, j, k, l, newsol, quad_evp = 0; int n, j, k, l, newsol, quad_evp = 0;
double tmp, d1, d2, abs, arg; double tmp, d1, d2, abs, arg;
complex_16 f, omega, omega2; complex_16 f, omega, omega2;
gMatrix *K = &DofData_P->M1; /* matrix associated with terms with no Dt nor Dt gMatrix *K =
Dt */ &DofData_P->M1; /* matrix associated with terms with no Dt nor DtDt */
gMatrix *L = &DofData_P->M2; /* matrix associated with Dt terms */ gMatrix *L = &DofData_P->M2; /* matrix associated with Dt terms */
gMatrix *M = &DofData_P->M3; /* matrix associated with DtDt terms */ gMatrix *M = &DofData_P->M3; /* matrix associated with DtDt terms */
gMatrix D; /* temp matrix for quadratic eigenvalue problem */ gMatrix D; /* temp matrix for quadratic eigenvalue problem */
/* Arpack parameters: see below for explanation */ /* Arpack parameters: see below for explanation */
int ido, nev, ncv, ldv, iparam[11], ipntr[14], lworkl, info, ldz; int ido, nev, ncv, ldv, iparam[11], ipntr[14], lworkl, info, ldz;
char bmat, *which, howmny; char bmat, *which, howmny;
double tol, *rwork; double tol, *rwork;
unsigned int rvec, *select; unsigned int rvec, *select;
complex_16 *resid, *v, *workd, *workl, *d, *z, sigma, *workev; complex_16 *resid, *v, *workd, *workl, *d, *z, sigma, *workev;
/* Warn if we are not in harmonic regime (we won't be able to compute/store /* Warn if we are not in harmonic regime (we won't be able to compute/store
complex eigenvectors) */ complex eigenvectors) */
if(Current.NbrHar != 2){ if(Current.NbrHar != 2) {
Message::Info("EigenSolve will only store the real part of the eigenvectors; Message::Info(
" "EigenSolve will only store the real part of the eigenvectors; "
"Define the system with \"Type Complex\" if this is an issue") "Define the system with \"Type Complex\" if this is an issue");
;
} }
#if defined(HAVE_PETSC) && !defined(PETSC_USE_COMPLEX) #if defined(HAVE_PETSC) && !defined(PETSC_USE_COMPLEX)
if(Current.NbrHar == 2){ if(Current.NbrHar == 2) {
Message::Warning("Using PETSc in real arithmetic for complex-simulated-real Message::Warning(
matrices"); "Using PETSc in real arithmetic for complex-simulated-real matrices");
} }
#endif #endif
/* Sanity checks */ /* Sanity checks */
if(DofData_P->Flag_Init[7] || DofData_P->Flag_Init[6] || if(DofData_P->Flag_Init[7] || DofData_P->Flag_Init[6] ||
DofData_P->Flag_Init[5] || DofData_P->Flag_Init[4]){ DofData_P->Flag_Init[5] || DofData_P->Flag_Init[4]) {
Message::Error("High order polynomial and non-linear EVP only available with Message::Error(
SLEPc"); "High order polynomial and non-linear EVP only available with SLEPc");
return; return;
} }
if(!DofData_P->Flag_Init[1] || !DofData_P->Flag_Init[3]){ if(!DofData_P->Flag_Init[1] || !DofData_P->Flag_Init[3]) {
Message::Error("No System available for EigenSolve: check 'DtDt' and 'Genera Message::Error("No System available for EigenSolve: check 'DtDt' and "
teSeparate'"); "'GenerateSeparate'");
return; return;
} }
/* Check if we have a "quadratic" evp (- w^2 M x + i w L x + K x = 0) */ /* Check if we have a "quadratic" evp (- w^2 M x + i w L x + K x = 0) */
if(DofData_P->Flag_Init[2]) if(DofData_P->Flag_Init[2]) quad_evp = 1;
quad_evp = 1;
/* Get eigenproblem parameters */ /* Get eigenproblem parameters */
EigenPar("eigen.par", &eigenpar); EigenPar("eigen.par", &eigenpar);
/* size of the system */ /* size of the system */
int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1; int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1;
n = DofData_P->NbrDof / incr; n = DofData_P->NbrDof / incr;
if(quad_evp) if(quad_evp) n *= 2;
n *= 2;
ido = 0; ido = 0;
/* Reverse communication flag. IDO must be zero on the first /* Reverse communication flag. IDO must be zero on the first
call to znaupd. IDO will be set internally to call to znaupd. IDO will be set internally to
indicate the type of operation to be performed. Control is indicate the type of operation to be performed. Control is
then given back to the calling routine which has the then given back to the calling routine which has the
responsibility to carry out the requested operation and call responsibility to carry out the requested operation and call
znaupd with the result. The operand is given in znaupd with the result. The operand is given in
WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)). WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)).
------------------------------------------------------------- -------------------------------------------------------------
skipping to change at line 317 skipping to change at line 323
After the initialization phase, when the routine is used in After the initialization phase, when the routine is used in
the "shift-and-invert" mode, the vector M * X is already the "shift-and-invert" mode, the vector M * X is already
available and does not need to be recomputed in forming OP*X. */ available and does not need to be recomputed in forming OP*X. */
bmat = 'I'; bmat = 'I';
/* BMAT specifies the type of the matrix B that defines the /* BMAT specifies the type of the matrix B that defines the
semi-inner product for the operator OP. semi-inner product for the operator OP.
BMAT = 'I' -> standard eigenvalue problem A*x = lambda*x BMAT = 'I' -> standard eigenvalue problem A*x = lambda*x
BMAT = 'G' -> generalized eigenvalue problem A*x = lambda*M*x */ BMAT = 'G' -> generalized eigenvalue problem A*x = lambda*M*x */
which = (char*)"LM"; which = (char *)"LM";
/* Which eigenvalues we want: /* Which eigenvalues we want:
SM = smallest magnitude ( magnitude = absolute value ) SM = smallest magnitude ( magnitude = absolute value )
LM = largest magnitude LM = largest magnitude
SR = smallest real part SR = smallest real part
LR = largest real part LR = largest real part
SI = smallest imaginary part SI = smallest imaginary part
LI = largest imaginary part */ LI = largest imaginary part */
nev = NumEigenvalues; nev = NumEigenvalues;
/* Number of eigenvalues of OP to be computed. 0 < NEV < N-1. /* Number of eigenvalues of OP to be computed. 0 < NEV < N-1.
Therefore, you'll be able to compute AT MOST n-2 eigenvalues! */ Therefore, you'll be able to compute AT MOST n-2 eigenvalues! */
/* sanity check */ /* sanity check */
if(nev >= n-1){ if(nev >= n - 1) {
Message::Warning("NumEigenvalues too large (%d < %d): setting to %d", nev, n Message::Warning("NumEigenvalues too large (%d < %d): setting to %d", nev,
-1, n-2); n - 1, n - 2);
nev = n-2; nev = n - 2;
} }
tol = eigenpar.prec; /* 1.e-4; */ tol = eigenpar.prec; /* 1.e-4; */
/* Stopping criteria: the relative accuracy of the Ritz value /* Stopping criteria: the relative accuracy of the Ritz value
is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I)) is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I))
where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex. where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex.
DEFAULT = dlamch('EPS') (machine precision as computed DEFAULT = dlamch('EPS') (machine precision as computed
by the LAPACK auxiliary subroutine dlamch). */ by the LAPACK auxiliary subroutine dlamch). */
resid = (complex_16*)Malloc(n * sizeof(complex_16)); resid = (complex_16 *)Malloc(n * sizeof(complex_16));
/* On INPUT: /* On INPUT:
If INFO .EQ. 0, a random initial residual vector is used. If INFO .EQ. 0, a random initial residual vector is used.
If INFO .NE. 0, RESID contains the initial residual vector, If INFO .NE. 0, RESID contains the initial residual vector,
possibly from a previous run. possibly from a previous run.
On OUTPUT: On OUTPUT:
RESID contains the final residual vector. */ RESID contains the final residual vector. */
ncv = eigenpar.size; /* Rule of thumb: NumEigenvalues * 2; */ ncv = eigenpar.size; /* Rule of thumb: NumEigenvalues * 2; */
/* Number of columns of the matrix V. NCV must satisfy the two /* Number of columns of the matrix V. NCV must satisfy the two
inequalities 1 <= NCV-NEV and NCV <= N. inequalities 1 <= NCV-NEV and NCV <= N.
This will indicate how many Arnoldi vectors are generated This will indicate how many Arnoldi vectors are generated
at each iteration. After the startup phase in which NEV at each iteration. After the startup phase in which NEV
Arnoldi vectors are generated, the algorithm generates Arnoldi vectors are generated, the algorithm generates
approximately NCV-NEV Arnoldi vectors at each subsequent update approximately NCV-NEV Arnoldi vectors at each subsequent update
iteration. Most of the cost in generating each Arnoldi vector is iteration. Most of the cost in generating each Arnoldi vector is
in the matrix-vector operation OP*x. */ in the matrix-vector operation OP*x. */
/* sanity checks */ /* sanity checks */
if(ncv <= nev){ if(ncv <= nev) {
Message::Warning("Krylov space size too small (%d <= %d), setting to %d", nc Message::Warning("Krylov space size too small (%d <= %d), setting to %d",
v, nev, nev*2); ncv, nev, nev * 2);
ncv = nev * 2; ncv = nev * 2;
} }
if(ncv > n){ if(ncv > n) {
Message::Warning("Krylov space size too large (%d > %d), setting to %d", ncv Message::Warning("Krylov space size too large (%d > %d), setting to %d",
, n, n); ncv, n, n);
ncv = n; ncv = n;
} }
v = (complex_16*)Malloc(n * ncv * sizeof(complex_16)); v = (complex_16 *)Malloc(n * ncv * sizeof(complex_16));
/* At the end of calculations, here will be stored the Arnoldi basis /* At the end of calculations, here will be stored the Arnoldi basis
vectors */ vectors */
ldv = n; ldv = n;
/* Leading dimension of "v". In our case, the number of lines of /* Leading dimension of "v". In our case, the number of lines of
"v". */ "v". */
iparam[0] = 1; iparam[0] = 1;
iparam[1] = 0; iparam[1] = 0;
iparam[2] = 10000; iparam[2] = 10000;
skipping to change at line 467 skipping to change at line 476
IPNTR(9): pointer to the NCV RITZ values of the IPNTR(9): pointer to the NCV RITZ values of the
original system. original system.
IPNTR(10): Not Used IPNTR(10): Not Used
IPNTR(11): pointer to the NCV corresponding error bounds. IPNTR(11): pointer to the NCV corresponding error bounds.
IPNTR(12): pointer to the NCV by NCV upper triangular IPNTR(12): pointer to the NCV by NCV upper triangular
Schur matrix for H. Schur matrix for H.
IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors
of the upper Hessenberg matrix H. Only referenced by of the upper Hessenberg matrix H. Only referenced by
zneupd if RVEC = .TRUE. See Remark 2 below. */ zneupd if RVEC = .TRUE. See Remark 2 below. */
workd = (complex_16*)Malloc(3 * n * sizeof(complex_16)); workd = (complex_16 *)Malloc(3 * n * sizeof(complex_16));
/* Distributed array to be used in the basic Arnoldi iteration /* Distributed array to be used in the basic Arnoldi iteration
for reverse communication. The user should not use WORKD for reverse communication. The user should not use WORKD
as temporary workspace during the iteration !!!!!!!!!! as temporary workspace during the iteration !!!!!!!!!!
See Data Distribution Note below. */ See Data Distribution Note below. */
lworkl = 3*ncv*ncv + 5*ncv; lworkl = 3 * ncv * ncv + 5 * ncv;
/* Dimension of the "workl" vector (see below). On must have: /* Dimension of the "workl" vector (see below). On must have:
lworkl >= 3*ncv*ncv + 5*ncv */ lworkl >= 3*ncv*ncv + 5*ncv */
workl = (complex_16*)Malloc(lworkl * sizeof(complex_16)); workl = (complex_16 *)Malloc(lworkl * sizeof(complex_16));
/* Private (replicated) array on each PE or array allocated on /* Private (replicated) array on each PE or array allocated on
the front end. See Data Distribution Note below. */ the front end. See Data Distribution Note below. */
rwork = (double*)Malloc(ncv * sizeof(double)); rwork = (double *)Malloc(ncv * sizeof(double));
/* Used as workspace */ /* Used as workspace */
info = 0; info = 0;
/* If INFO .EQ. 0, a randomly initial residual vector is used. /* If INFO .EQ. 0, a randomly initial residual vector is used.
If INFO .NE. 0, RESID contains the initial residual vector, If INFO .NE. 0, RESID contains the initial residual vector,
possibly from a previous run. possibly from a previous run.
Error flag on output. Error flag on output.
= 0: Normal exit. = 0: Normal exit.
= 1: Maximum number of iterations taken. = 1: Maximum number of iterations taken.
All possible eigenvalues of OP has been found. IPARAM(5) All possible eigenvalues of OP has been found. IPARAM(5)
skipping to change at line 525 skipping to change at line 534
IPARAM(5) returns the size of the current Arnoldi IPARAM(5) returns the size of the current Arnoldi
factorization. */ factorization. */
rvec = 1; /* .true. */ rvec = 1; /* .true. */
/* If we want Ritz vectors to be computed as well. */ /* If we want Ritz vectors to be computed as well. */
howmny = 'A'; howmny = 'A';
/* What do we want: Ritz or Schur vectors? For Schur, choose: howmny /* What do we want: Ritz or Schur vectors? For Schur, choose: howmny
= 'P' */ = 'P' */
select = (unsigned int*)Malloc(ncv * sizeof(unsigned int)); select = (unsigned int *)Malloc(ncv * sizeof(unsigned int));
/* Internal workspace */ /* Internal workspace */
d = (complex_16*)Malloc(nev * sizeof(complex_16)); d = (complex_16 *)Malloc(nev * sizeof(complex_16));
/* Vector containing the "nev" eigenvalues computed. /* Vector containing the "nev" eigenvalues computed.
VERY IMPORTANT: on line 69 of zneupd.f they say it should be nev+1; VERY IMPORTANT: on line 69 of zneupd.f they say it should be nev+1;
this is wrong, for see line 283 where it is declared as d(nev) */ this is wrong, for see line 283 where it is declared as d(nev) */
z = (complex_16*)Malloc(n * nev * sizeof(complex_16)); z = (complex_16 *)Malloc(n * nev * sizeof(complex_16));
/* On exit, if RVEC = .TRUE. and HOWMNY = 'A', then the columns of /* On exit, if RVEC = .TRUE. and HOWMNY = 'A', then the columns of
Z represents approximate eigenvectors (Ritz vectors) corresponding Z represents approximate eigenvectors (Ritz vectors) corresponding
to the NCONV=IPARAM(5) Ritz values for eigensystem to the NCONV=IPARAM(5) Ritz values for eigensystem
A*z = lambda*B*z. A*z = lambda*B*z.
If RVEC = .FALSE. or HOWMNY = 'P', then Z is NOT REFERENCED. If RVEC = .FALSE. or HOWMNY = 'P', then Z is NOT REFERENCED.
NOTE: If if RVEC = .TRUE. and a Schur basis is not required, NOTE: If if RVEC = .TRUE. and a Schur basis is not required,
the array Z may be set equal to first NEV+1 columns of the Arnoldi the array Z may be set equal to first NEV+1 columns of the Arnoldi
basis array V computed by ZNAUPD. In this case the Arnoldi basis basis array V computed by ZNAUPD. In this case the Arnoldi basis
will be destroyed and overwritten with the eigenvector basis. */ will be destroyed and overwritten with the eigenvector basis. */
ldz = n; ldz = n;
/* Leading dimension of "z". In our case, the number of lines of "z". */ /* Leading dimension of "z". In our case, the number of lines of "z". */
sigma.re = 0.; sigma.re = 0.;
sigma.im = 0.; sigma.im = 0.;
/* The shift. Not used in this case: we deal with the shift "by /* The shift. Not used in this case: we deal with the shift "by
hand". */ hand". */
workev = (complex_16*)Malloc(2 * ncv * sizeof(complex_16)); workev = (complex_16 *)Malloc(2 * ncv * sizeof(complex_16));
/* Workspace */ /* Workspace */
if(bmat != 'I' || iparam[6] != 1){ if(bmat != 'I' || iparam[6] != 1) {
Message::Error("General and/or shift-invert mode should not be used"); Message::Error("General and/or shift-invert mode should not be used");
return; return;
} }
/* Create temp vectors and matrices and apply shift. Warning: with /* Create temp vectors and matrices and apply shift. Warning: with
PETSc, the shifting can be very slow if the masks are very PETSc, the shifting can be very slow if the masks are very
different, for example if we are in real arithmetic and have one different, for example if we are in real arithmetic and have one
real matrix and one complex "simulated-real" matrix */ real matrix and one complex "simulated-real" matrix */
if(!quad_evp){ if(!quad_evp) {
LinAlg_CreateVector(&v1, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&v1, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CreateVector(&v2, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&v2, &DofData_P->Solver, DofData_P->NbrDof);
/* K = K - shift * M */ /* K = K - shift * M */
LinAlg_AddMatrixProdMatrixDouble(K, M, -shift_r, K) ; LinAlg_AddMatrixProdMatrixDouble(K, M, -shift_r, K);
} }
else{ else {
/* This is an explanation of our approach to a quadratic /* This is an explanation of our approach to a quadratic
eigenvalue problem i.e. - w^2 M x + i w L x + K x = 0. This eigenvalue problem i.e. - w^2 M x + i w L x + K x = 0. This
system is equivalent to (y = i w x) and (i w M y + i w L x + K system is equivalent to (y = i w x) and (i w M y + i w L x + K
x = 0), or, in matrix form: x = 0), or, in matrix form:
| L M | |x| |-K 0 | |x| | L M | |x| |-K 0 | |x|
| I 0 | |y| iw = | 0 I | |y| , or | I 0 | |y| iw = | 0 I | |y| , or
|x| |x| |x| |x|
A |y| iw = B |y|. A |y| iw = B |y|.
skipping to change at line 602 skipping to change at line 611
matrices but the Arpack vector is of size n=2*N. Only the x matrices but the Arpack vector is of size n=2*N. Only the x
part of the (x,y) eigenvectors are retained as physical part of the (x,y) eigenvectors are retained as physical
solutions. */ solutions. */
LinAlg_CreateVector(&x, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&x, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CreateVector(&y, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&y, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CreateVector(&v1, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&v1, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CreateVector(&w1, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&w1, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CreateVector(&w2, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&w2, &DofData_P->Solver, DofData_P->NbrDof);
LinAlg_CreateMatrix(&D, &DofData_P->Solver, DofData_P->NbrDof, DofData_P->Nb LinAlg_CreateMatrix(&D, &DofData_P->Solver, DofData_P->NbrDof,
rDof); DofData_P->NbrDof);
/* D = -(shift^2 * M + shift * L + K) */ /* D = -(shift^2 * M + shift * L + K) */
LinAlg_CopyMatrix(M, &D); LinAlg_CopyMatrix(M, &D);
LinAlg_AddMatrixProdMatrixDouble(L, &D, shift_r, &D); LinAlg_AddMatrixProdMatrixDouble(L, &D, shift_r, &D);
LinAlg_AddMatrixProdMatrixDouble(K, &D, shift_r, &D); LinAlg_AddMatrixProdMatrixDouble(K, &D, shift_r, &D);
LinAlg_ProdMatrixDouble(&D, -1., &D); LinAlg_ProdMatrixDouble(&D, -1., &D);
} }
/* Keep calling znaupd again and again until ido == 99 */ /* Keep calling znaupd again and again until ido == 99 */
k = 0; k = 0;
do { do {
znaupd_(&ido, &bmat, &n, which, &nev, &tol, resid, &ncv, v, &ldv, iparam, znaupd_(&ido, &bmat, &n, which, &nev, &tol, resid, &ncv, v, &ldv, iparam,
ipntr, workd, workl, &lworkl, rwork, &info); ipntr, workd, workl, &lworkl, rwork, &info);
if(ido == 1 || ido == -1){ if(ido == 1 || ido == -1) {
Message::Info("Arpack iteration %d", k+1); Message::Info("Arpack iteration %d", k + 1);
if(!quad_evp){ if(!quad_evp) {
Arpack2GetDP(n, &workd[ipntr[0]-1], &v1); Arpack2GetDP(n, &workd[ipntr[0] - 1], &v1);
LinAlg_ProdMatrixVector(M, &v1, &v2); LinAlg_ProdMatrixVector(M, &v1, &v2);
if(!k) if(!k)
LinAlg_Solve(K, &v2, &DofData_P->Solver, &v1); LinAlg_Solve(K, &v2, &DofData_P->Solver, &v1);
else else
LinAlg_SolveAgain(K, &v2, &DofData_P->Solver, &v1); LinAlg_SolveAgain(K, &v2, &DofData_P->Solver, &v1);
GetDP2Arpack(&v1, &workd[ipntr[1]-1]); GetDP2Arpack(&v1, &workd[ipntr[1] - 1]);
} }
else{ else {
Arpack2GetDPSplit(n, &workd[ipntr[0]-1], &x, &y); Arpack2GetDPSplit(n, &workd[ipntr[0] - 1], &x, &y);
LinAlg_ProdMatrixVector(M, &y, &w2); LinAlg_ProdMatrixVector(M, &y, &w2);
LinAlg_ProdMatrixVector(L, &x, &v1); LinAlg_ProdMatrixVector(L, &x, &v1);
LinAlg_AddVectorVector(&v1, &w2, &v1); LinAlg_AddVectorVector(&v1, &w2, &v1);
LinAlg_ProdMatrixVector(M, &x, &w1); LinAlg_ProdMatrixVector(M, &x, &w1);
LinAlg_AddVectorProdVectorDouble(&v1, &w1, shift_r, &v1); LinAlg_AddVectorProdVectorDouble(&v1, &w1, shift_r, &v1);
if(!k) if(!k)
LinAlg_Solve(&D, &v1, &DofData_P->Solver, &w1); LinAlg_Solve(&D, &v1, &DofData_P->Solver, &w1);
else else
LinAlg_SolveAgain(&D, &v1, &DofData_P->Solver, &w1); LinAlg_SolveAgain(&D, &v1, &DofData_P->Solver, &w1);
LinAlg_ProdMatrixVector(K, &x, &v1); LinAlg_ProdMatrixVector(K, &x, &v1);
LinAlg_ProdVectorDouble(&v1, -1., &v1); LinAlg_ProdVectorDouble(&v1, -1., &v1);
LinAlg_AddVectorProdVectorDouble(&v1, &w2, shift_r, &v1); LinAlg_AddVectorProdVectorDouble(&v1, &w2, shift_r, &v1);
LinAlg_SolveAgain(&D, &v1, &DofData_P->Solver, &w2); LinAlg_SolveAgain(&D, &v1, &DofData_P->Solver, &w2);
GetDP2ArpackMerge(&w1, &w2, &workd[ipntr[1]-1]); GetDP2ArpackMerge(&w1, &w2, &workd[ipntr[1] - 1]);
} }
k++; k++;
} }
else if(ido == 99){ else if(ido == 99) {
/* We're done! */ /* We're done! */
break; break;
} }
else{ else {
Message::Info("Arpack code = %d (ignored)", info); Message::Info("Arpack code = %d (ignored)", info);
} }
} while (1); } while(1);
Message::Info("Arpack required %d iterations", k); Message::Info("Arpack required %d iterations", k);
/* Testing for errors */ /* Testing for errors */
if(info == 0){ if(info == 0) { /* OK */
/* OK */
} }
else if(info == 1){ else if(info == 1) {
Message::Warning("Maxmimum number of iteration reached in EigenSolve"); Message::Warning("Maxmimum number of iteration reached in EigenSolve");
} }
else if(info == 2){ else if(info == 2) {
Message::Warning("No shifts could be applied during a cycle of the"); Message::Warning("No shifts could be applied during a cycle of the");
Message::Warning("Implicitly restarted Arnoldi iteration. One possibility"); Message::Warning("Implicitly restarted Arnoldi iteration. One possibility");
Message::Warning("is to increase the size of NCV relative to NEV."); Message::Warning("is to increase the size of NCV relative to NEV.");
} }
else if(info < 0){ else if(info < 0) {
Message::Error("Arpack code = %d", info); Message::Error("Arpack code = %d", info);
} }
else{ else {
Message::Warning("Arpack code = %d (unknown)", info); Message::Warning("Arpack code = %d (unknown)", info);
} }
/* Call to zneupd for post-processing */ /* Call to zneupd for post-processing */
zneupd_(&rvec, &howmny, select, d, z, &ldz, &sigma, workev, &bmat, &n, which, zneupd_(&rvec, &howmny, select, d, z, &ldz, &sigma, workev, &bmat, &n, which,
&nev, &tol, resid, &ncv, v, &ldv, iparam, ipntr, workd, workl, &lworkl, &nev, &tol, resid, &ncv, v, &ldv, iparam, ipntr, workd, workl,
rwork, &info); &lworkl, rwork, &info);
/* Test for errors */ /* Test for errors */
if(info != 0) if(info != 0)
Message::Error("Arpack code = %d (eigenvector post-processing)", info); Message::Error("Arpack code = %d (eigenvector post-processing)", info);
/* Compute the unshifted eigenvalues and print them, and store the /* Compute the unshifted eigenvalues and print them, and store the
associated eigenvectors */ associated eigenvectors */
newsol = 0; newsol = 0;
for (k = 0; k < nev; k++){ for(k = 0; k < nev; k++) {
/* Unshift the eigenvalues */ /* Unshift the eigenvalues */
tmp = SQU(d[k].re) + SQU(d[k].im); tmp = SQU(d[k].re) + SQU(d[k].im);
d[k].re = shift_r + d[k].re/tmp; d[k].re = shift_r + d[k].re / tmp;
d[k].im = shift_i - d[k].im/tmp; d[k].im = shift_i - d[k].im / tmp;
if(!quad_evp){ if(!quad_evp) {
/* Eigenvalue = omega^2 */ /* Eigenvalue = omega^2 */
omega2.re = d[k].re; omega2.re = d[k].re;
omega2.im = d[k].im; omega2.im = d[k].im;
abs = sqrt(SQU(omega2.re) + SQU(omega2.im)); abs = sqrt(SQU(omega2.re) + SQU(omega2.im));
arg = atan2(omega2.im, omega2.re); arg = atan2(omega2.im, omega2.re);
omega.re = sqrt(abs) * cos(0.5*arg); omega.re = sqrt(abs) * cos(0.5 * arg);
omega.im = sqrt(abs) * sin(0.5*arg); omega.im = sqrt(abs) * sin(0.5 * arg);
f.re = omega.re / TWO_PI; f.re = omega.re / TWO_PI;
f.im = omega.im / TWO_PI; f.im = omega.im / TWO_PI;
} }
else{ else {
/* Eigenvalue = i*omega */ /* Eigenvalue = i*omega */
omega.re = d[k].im; omega.re = d[k].im;
omega.im = -d[k].re; omega.im = -d[k].re;
omega2.re = SQU(omega.re) - SQU(omega.im); omega2.re = SQU(omega.re) - SQU(omega.im);
omega2.im = 2. * omega.re * omega.im; omega2.im = 2. * omega.re * omega.im;
f.re = omega.re / TWO_PI; f.re = omega.re / TWO_PI;
f.im = omega.im / TWO_PI; f.im = omega.im / TWO_PI;
} }
Message::Info("Eigenvalue %03d: w^2 = %.12e %s %.12e * i", Message::Info("Eigenvalue %03d: w^2 = %.12e %s %.12e * i", k + 1, omega2.re,
k+1, omega2.re, (omega2.im > 0) ? "+" : "-", (omega2.im > 0) ? "+" : "-",
(omega2.im > 0) ? omega2.im : -omega2.im); (omega2.im > 0) ? omega2.im : -omega2.im);
Message::Info(" w = %.12e %s %.12e * i", Message::Info(" w = %.12e %s %.12e * i", omega.re,
omega.re, (omega.im > 0) ? "+" : "-", (omega.im > 0) ? "+" : "-",
(omega.im > 0) ? omega.im : -omega.im); (omega.im > 0) ? omega.im : -omega.im);
Message::Info(" f = %.12e %s %.12e * i", Message::Info(" f = %.12e %s %.12e * i", f.re,
f.re, (f.im > 0) ? "+" : "-", (f.im > 0) ? f.im : -f.im); (f.im > 0) ? "+" : "-", (f.im > 0) ? f.im : -f.im);
/* Update the current value of Time and TimeImag so that /* Update the current value of Time and TimeImag so that
$EigenvalueReal and $EigenvalueImag are up-to-date */ $EigenvalueReal and $EigenvalueImag are up-to-date */
Current.Time = omega.re; Current.Time = omega.re;
Current.TimeImag = omega.im; Current.TimeImag = omega.im;
// test filter expression and continue without storing if false // test filter expression and continue without storing if false
if(FilterExpressionIndex >= 0){ if(FilterExpressionIndex >= 0) {
struct Value val; struct Value val;
Get_ValueOfExpressionByIndex(FilterExpressionIndex, NULL, 0., 0., 0., &val Get_ValueOfExpressionByIndex(FilterExpressionIndex, NULL, 0., 0., 0.,
); &val);
if(!val.Val[0]){ if(!val.Val[0]) {
Message::Debug("Skipping eigenvalue %g + i * %g", omega.re, omega.im); Message::Debug("Skipping eigenvalue %g + i * %g", omega.re, omega.im);
continue; continue;
} }
} }
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Re(Omega)" Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
, "/Re(Omega)",
std::vector<double>(1, omega.re)); std::vector<double>(1, omega.re));
Message::AddOnelabNumberChoice(Message::GetOnelabClientName() + "/Im(Omega)" Message::AddOnelabNumberChoice(Message::GetOnelabClientName() +
, "/Im(Omega)",
std::vector<double>(1, omega.im)); std::vector<double>(1, omega.im));
if(newsol) { if(newsol) {
/* Create new solution */ /* Create new solution */
Solution_S.TimeFunctionValues = NULL; Solution_S.TimeFunctionValues = NULL;
LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof); LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof);
List_Add(DofData_P->Solutions, &Solution_S); List_Add(DofData_P->Solutions, &Solution_S);
DofData_P->CurrentSolution = (struct Solution*) DofData_P->CurrentSolution = (struct Solution *)List_Pointer(
List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1); DofData_P->Solutions, List_Nbr(DofData_P->Solutions) - 1);
} }
newsol = 1; newsol = 1;
DofData_P->CurrentSolution->Time = omega.re; DofData_P->CurrentSolution->Time = omega.re;
DofData_P->CurrentSolution->TimeImag = omega.im; DofData_P->CurrentSolution->TimeImag = omega.im;
DofData_P->CurrentSolution->TimeStep = (int)Current.TimeStep; DofData_P->CurrentSolution->TimeStep = (int)Current.TimeStep;
Free(DofData_P->CurrentSolution->TimeFunctionValues); Free(DofData_P->CurrentSolution->TimeFunctionValues);
DofData_P->CurrentSolution->TimeFunctionValues = NULL; DofData_P->CurrentSolution->TimeFunctionValues = NULL;
DofData_P->CurrentSolution->SolutionExist = 1; DofData_P->CurrentSolution->SolutionExist = 1;
int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1; int incr = (Current.NbrHar == 2) ? gCOMPLEX_INCREMENT : 1;
for(l = 0; l < DofData_P->NbrDof; l += incr){ for(l = 0; l < DofData_P->NbrDof; l += incr) {
j = l / incr; j = l / incr;
if(Current.NbrHar == 2){ if(Current.NbrHar == 2) {
LinAlg_SetComplexInVector(z[k*n+j].re, z[k*n+j].im, LinAlg_SetComplexInVector(z[k * n + j].re, z[k * n + j].im,
&DofData_P->CurrentSolution->x, l, l+1); &DofData_P->CurrentSolution->x, l, l + 1);
} }
else{ else {
LinAlg_SetDoubleInVector(z[k*n+j].re, LinAlg_SetDoubleInVector(z[k * n + j].re,
&DofData_P->CurrentSolution->x, l); &DofData_P->CurrentSolution->x, l);
} }
} }
LinAlg_AssembleVector(&DofData_P->CurrentSolution->x); LinAlg_AssembleVector(&DofData_P->CurrentSolution->x);
/* Arpack returns eigenvectors normalized in L-2 norm. Renormalize /* Arpack returns eigenvectors normalized in L-2 norm. Renormalize
them in L-infty norm so that the absolute value of the largest them in L-infty norm so that the absolute value of the largest
element is 1 */ element is 1 */
tmp = 0.; tmp = 0.;
for(l = 0; l < DofData_P->NbrDof; l += incr){ for(l = 0; l < DofData_P->NbrDof; l += incr) {
if(Current.NbrHar == 2){ if(Current.NbrHar == 2) {
LinAlg_GetComplexInVector(&d1, &d2, LinAlg_GetComplexInVector(&d1, &d2, &DofData_P->CurrentSolution->x, l,
&DofData_P->CurrentSolution->x, l, l+1); l + 1);
abs = sqrt(SQU(d1) + SQU(d2)); abs = sqrt(SQU(d1) + SQU(d2));
} }
else{ else {
LinAlg_GetDoubleInVector(&d1, LinAlg_GetDoubleInVector(&d1, &DofData_P->CurrentSolution->x, l);
&DofData_P->CurrentSolution->x, l);
abs = sqrt(SQU(d1)); abs = sqrt(SQU(d1));
} }
if(abs > tmp) tmp = abs; if(abs > tmp) tmp = abs;
} }
if(tmp > 1.e-16) if(tmp > 1.e-16)
LinAlg_ProdVectorDouble(&DofData_P->CurrentSolution->x, 1./tmp, LinAlg_ProdVectorDouble(&DofData_P->CurrentSolution->x, 1. / tmp,
&DofData_P->CurrentSolution->x); &DofData_P->CurrentSolution->x);
/* Increment the global timestep counter so that a future /* Increment the global timestep counter so that a future
GenerateSystem knows which solutions exist */ GenerateSystem knows which solutions exist */
Current.TimeStep += 1.; Current.TimeStep += 1.;
} }
/* Deallocate */ /* Deallocate */
if(!quad_evp){ if(!quad_evp) {
LinAlg_DestroyVector(&v1); LinAlg_DestroyVector(&v1);
LinAlg_DestroyVector(&v2); LinAlg_DestroyVector(&v2);
} }
else{ else {
LinAlg_DestroyVector(&x); LinAlg_DestroyVector(&x);
LinAlg_DestroyVector(&y); LinAlg_DestroyVector(&y);
LinAlg_DestroyVector(&v1); LinAlg_DestroyVector(&v1);
LinAlg_DestroyVector(&w1); LinAlg_DestroyVector(&w1);
LinAlg_DestroyVector(&w2); LinAlg_DestroyVector(&w2);
LinAlg_DestroyMatrix(&D); LinAlg_DestroyMatrix(&D);
} }
Free(resid); Free(resid);
Free(v); Free(v);
 End of changes. 86 change blocks. 
197 lines changed or deleted 189 lines changed or added
To the C++ Programs section of the getdp source changes report or the top of this page
Mainly generated by pkgdiff using rfcdiff
Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)