SALSA Analysis Modules
|
Dense routines from Lapack: eigenvalue and schur stuff. More...
#include <stdlib.h>
#include "anamod.h"
#include "anamodsalsamodules.h"
#include "petscerror.h"
#include "petscmat.h"
Go to the source code of this file.
Defines | |
#define | ABS(x) (x>0?x:-x) |
#define | EIGENVALUE(re, im) |
Functions | |
void | LAPACK_DGEESX (const char *jobvs, const char *sort, char *sel, const char *sense, const int *n, double *a, const int *lda, int *sdim, double *wr, double *wi, double *vs, const int *ldvs, double *rconde, double *rcondv, double *work, const int *lwork, int *iwork, const int *liwork, int *bwork, int *info) |
static PetscErrorCode | eigenvaluecomp (Mat A) |
static PetscErrorCode | Departure (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg) |
static PetscErrorCode | MaxEVbyMagRe (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg) |
static PetscErrorCode | MaxEVbyMagIm (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg) |
static PetscErrorCode | MinEVbyMagRe (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg) |
static PetscErrorCode | MinEVbyMagIm (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg) |
static PetscErrorCode | MaxEVbyRealRe (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg) |
static PetscErrorCode | MaxEVbyRealIm (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg) |
static PetscErrorCode | MaxEVbyImRe (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg) |
static PetscErrorCode | MaxEVbyImIm (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscBool *flg) |
PetscErrorCode | RegisterLapackModules () |
Dense routines from Lapack: eigenvalue and schur stuff.
Activate this module with
PetscErrorCode RegisterLapackModules();
Compute these elements with
ComputeQuantity("lapack",<element>,A,(AnalysisItem*)&res,&flg);
Available elements are:
Definition in file lapack.c.
#define ABS | ( | x | ) | (x>0?x:-x) |
Referenced by eigenvaluecomp().
#define EIGENVALUE | ( | re, | |
im | |||
) |
{ \ PetscReal mag = sqrt(re*re+im*im); \ if (!maxbyreset) {maxbyreset = 1; maxbyrere = re; maxbyreim = im; } \ else if (re>maxbyrere) { maxbyrere = re; maxbyreim = im; } \ if (!maxbyimset) {maxbyimset = 1; maxbyimre = re; maxbyimim = im; } \ else if (im>maxbyimre) { maxbyimre = re; maxbyimim = im; } \ if (!maxbymagset) {maxbymagset = 1; maxmag = mag; \ maxbymagre = re; maxbymagim = im; } \ else if (mag>maxmag) { maxbymagre = re; maxbymagim = im; } \ if (!minbymagset) {minbymagset = 1; minmag = mag; \ minbymagre = re; minbymagim = im; } \ else if (mag<minmag) { minbymagre = re; minbymagim = im; } \ }
Referenced by eigenvaluecomp().
static PetscErrorCode Departure | ( | AnaModNumericalProblem | prob, |
AnalysisItem * | rv, | ||
int * | lv, | ||
PetscBool * | flg | ||
) | [static] |
Definition at line 196 of file lapack.c.
References eigenvaluecomp(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.
Referenced by RegisterLapackModules().
{ Mat A = (Mat)prob; PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; PetscFunctionBegin; ierr = GetDataID("lapack","departure",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else { ierr = eigenvaluecomp(A); CHKERRQ(ierr); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed departure by now"); } PetscFunctionReturn(0); }
static PetscErrorCode eigenvaluecomp | ( | Mat | A | ) | [static] |
This is the computational routine for lapack eigenvalue calculation.
Definition at line 41 of file lapack.c.
References ABS, EIGENVALUE, GetDataID(), HASTOEXIST, id, and LAPACK_DGEESX().
Referenced by Departure(), MaxEVbyImIm(), MaxEVbyImRe(), MaxEVbyMagIm(), MaxEVbyMagRe(), MaxEVbyRealIm(), MaxEVbyRealRe(), MinEVbyMagIm(), and MinEVbyMagRe().
{ MPI_Comm comm; Mat Adense; int maxbyreset=0,maxbyimset=0,maxbymagset=0,minbymagset=0; PetscReal dep=0, maxbymagre,maxbymagim, minbymagre,minbymagim, maxbyrere,maxbyreim, maxbyimre,maxbyimim, maxmag,minmag; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); ierr = MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&Adense); CHKERRQ(ierr); { char JOBVS = 'N', SORT = 'N', SENSE = 'N', SELECT='X'; int N,LDA,SDIM=0,LDVS=1,LWORK,LIWORK,*IWORK,INFO; PetscScalar *AA,*VS=NULL,*RCONDE=NULL,*RCONDV=NULL; PetscReal *WR,*WI,*WORK,*BWORK=NULL; ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr); LDA = N; ierr = MatGetArray(Adense,&AA); CHKERRQ(ierr); ierr = PetscMalloc(N*sizeof(PetscReal),&WR); CHKERRQ(ierr); ierr = PetscMalloc(N*sizeof(PetscReal),&WI); CHKERRQ(ierr); LWORK = 3*N+1; ierr = PetscMalloc(LWORK*sizeof(PetscReal),&WORK); CHKERRQ(ierr); LIWORK = 10; ierr = PetscMalloc(LIWORK*sizeof(int),&IWORK); CHKERRQ(ierr); CHKMEMQ; LAPACK_DGEESX(&JOBVS, &SORT, &SELECT, &SENSE, &N, AA, &LDA, &SDIM, WR, WI, VS, &LDVS, RCONDE, RCONDV, WORK, &LWORK, IWORK, &LIWORK, (int*)BWORK, &INFO); if (INFO<0) { INFO = -INFO; SETERRQ1(MPI_COMM_WORLD,1,"DGEESX Problem with parameter %d",INFO); } else if (INFO>0) { SETERRQ1(MPI_COMM_WORLD,1,"DGEESX error %d",INFO); } CHKMEMQ; { int i,j,rind=0; PetscReal f=0; for (i=0; i<N; i++) { #define ABS(x) (x>0?x:-x) #define EIGENVALUE(re,im) { \ PetscReal mag = sqrt(re*re+im*im); \ if (!maxbyreset) {maxbyreset = 1; maxbyrere = re; maxbyreim = im; } \ else if (re>maxbyrere) { maxbyrere = re; maxbyreim = im; } \ if (!maxbyimset) {maxbyimset = 1; maxbyimre = re; maxbyimim = im; } \ else if (im>maxbyimre) { maxbyimre = re; maxbyimim = im; } \ if (!maxbymagset) {maxbymagset = 1; maxmag = mag; \ maxbymagre = re; maxbymagim = im; } \ else if (mag>maxmag) { maxbymagre = re; maxbymagim = im; } \ if (!minbymagset) {minbymagset = 1; minmag = mag; \ minbymagre = re; minbymagim = im; } \ else if (mag<minmag) { minbymagre = re; minbymagim = im; } \ } PetscReal dd = AA[rind],db = ABS(AA[rind+1]),dt; dt = 1.e-14 * ABS(dd); /*printf("row %d, diagonal %e, below %e, test %e\n",i,dd,db,dt);*/ if (db<dt) { /* case of 1x1 eigenvalue block */ int ind=rind+N; for (j=i+1; j<N; j++) { PetscReal aa = AA[ind]; f += aa*aa; ind += N; } EIGENVALUE(dd,0.); } else { /* case of 2x2 eigenvalue block */ int ind; PetscReal dd = AA[rind],do1=AA[rind+1],do2=AA[rind+N],dn=AA[rind+N+1]; if (ABS((dd-dn)/dd)>1.e-13) SETERRQ4(PETSC_COMM_SELF,1,"Strange 2x2 block %e,%e,%e,%e\n",dd,do1,do2,dn); EIGENVALUE(dd,sqrt(do1*do2)); ind = rind+2*N; for (j=i+2; j<N; j++) { PetscReal aa = AA[ind]; f += aa*aa; ind += N; } i++; rind += N+1; ind = rind+N; for (j=i+1; j<N; j++) { PetscReal aa = AA[ind]; f += aa*aa; ind += N; } } rind += N+1; } dep = sqrt(f); } ierr = PetscFree(WR); CHKERRQ(ierr); ierr = PetscFree(WI); CHKERRQ(ierr); } ierr = MatDestroy(&Adense); CHKERRQ(ierr); { int id; PetscBool has; ierr = GetDataID("lapack","departure",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataSetScalar ((PetscObject)A,id,dep); CHKERRQ(ierr); ierr = GetDataID ("lapack","lambda-max-by-magnitude-re",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataSetReal ((PetscObject)A,id,maxbymagre); CHKERRQ(ierr); ierr = GetDataID ("lapack","lambda-max-by-magnitude-im",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataSetReal ((PetscObject)A,id,maxbymagim); CHKERRQ(ierr); ierr = GetDataID ("lapack","lambda-min-by-magnitude-re",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataSetReal ((PetscObject)A,id,minbymagre); CHKERRQ(ierr); ierr = GetDataID ("lapack","lambda-min-by-magnitude-im",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataSetReal ((PetscObject)A,id,minbymagim); CHKERRQ(ierr); ierr = GetDataID ("lapack","lambda-max-by-real-part-re",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataSetReal ((PetscObject)A,id,maxbyrere); CHKERRQ(ierr); ierr = GetDataID ("lapack","lambda-max-by-real-part-im",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataSetReal ((PetscObject)A,id,maxbyreim); CHKERRQ(ierr); ierr = GetDataID ("lapack","lambda-max-by-im-part-re",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataSetReal ((PetscObject)A,id,maxbyimre); CHKERRQ(ierr); ierr = GetDataID ("lapack","lambda-max-by-im-part-im",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataSetReal ((PetscObject)A,id,maxbyimim); CHKERRQ(ierr); } PetscFunctionReturn(0); }
void LAPACK_DGEESX | ( | const char * | jobvs, |
const char * | sort, | ||
char * | sel, | ||
const char * | sense, | ||
const int * | n, | ||
double * | a, | ||
const int * | lda, | ||
int * | sdim, | ||
double * | wr, | ||
double * | wi, | ||
double * | vs, | ||
const int * | ldvs, | ||
double * | rconde, | ||
double * | rcondv, | ||
double * | work, | ||
const int * | lwork, | ||
int * | iwork, | ||
const int * | liwork, | ||
int * | bwork, | ||
int * | info | ||
) |
Referenced by eigenvaluecomp().
static PetscErrorCode MaxEVbyImIm | ( | AnaModNumericalProblem | prob, |
AnalysisItem * | rv, | ||
int * | lv, | ||
PetscBool * | flg | ||
) | [static] |
Definition at line 397 of file lapack.c.
References eigenvaluecomp(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.
Referenced by RegisterLapackModules().
{ Mat A = (Mat)prob; PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; PetscFunctionBegin; ierr = GetDataID ("lapack","lambda-max-by-im-part-im",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataGetReal ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) rv->r = v; else { ierr = eigenvaluecomp(A); CHKERRQ(ierr); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-im-part-im by now"); } PetscFunctionReturn(0); }
static PetscErrorCode MaxEVbyImRe | ( | AnaModNumericalProblem | prob, |
AnalysisItem * | rv, | ||
int * | lv, | ||
PetscBool * | flg | ||
) | [static] |
Definition at line 372 of file lapack.c.
References eigenvaluecomp(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.
Referenced by RegisterLapackModules().
{ Mat A = (Mat)prob; PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; PetscFunctionBegin; ierr = GetDataID ("lapack","lambda-max-by-im-part-re",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataGetReal ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) rv->r = v; else { ierr = eigenvaluecomp(A); CHKERRQ(ierr); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-im-part-re by now"); } PetscFunctionReturn(0); }
static PetscErrorCode MaxEVbyMagIm | ( | AnaModNumericalProblem | prob, |
AnalysisItem * | rv, | ||
int * | lv, | ||
PetscBool * | flg | ||
) | [static] |
Definition at line 247 of file lapack.c.
References eigenvaluecomp(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.
Referenced by RegisterLapackModules().
{ Mat A = (Mat)prob; PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; PetscFunctionBegin; ierr = GetDataID ("lapack","lambda-max-by-magnitude-im",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataGetReal ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) rv->r = v; else { ierr = eigenvaluecomp(A); CHKERRQ(ierr); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-magnitude-im by now"); } PetscFunctionReturn(0); }
static PetscErrorCode MaxEVbyMagRe | ( | AnaModNumericalProblem | prob, |
AnalysisItem * | rv, | ||
int * | lv, | ||
PetscBool * | flg | ||
) | [static] |
Definition at line 222 of file lapack.c.
References eigenvaluecomp(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.
Referenced by RegisterLapackModules().
{ Mat A = (Mat)prob; PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; PetscFunctionBegin; ierr = GetDataID ("lapack","lambda-max-by-magnitude-re",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataGetReal ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) rv->r = v; else { ierr = eigenvaluecomp(A); CHKERRQ(ierr); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-magnitude-re by now"); } PetscFunctionReturn(0); }
static PetscErrorCode MaxEVbyRealIm | ( | AnaModNumericalProblem | prob, |
AnalysisItem * | rv, | ||
int * | lv, | ||
PetscBool * | flg | ||
) | [static] |
Definition at line 347 of file lapack.c.
References eigenvaluecomp(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.
Referenced by RegisterLapackModules().
{ Mat A = (Mat)prob; PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; PetscFunctionBegin; ierr = GetDataID ("lapack","lambda-max-by-real-part-im",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataGetReal ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) rv->r = v; else { ierr = eigenvaluecomp(A); CHKERRQ(ierr); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-real-part-im by now"); } PetscFunctionReturn(0); }
static PetscErrorCode MaxEVbyRealRe | ( | AnaModNumericalProblem | prob, |
AnalysisItem * | rv, | ||
int * | lv, | ||
PetscBool * | flg | ||
) | [static] |
Definition at line 322 of file lapack.c.
References eigenvaluecomp(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.
Referenced by RegisterLapackModules().
{ Mat A = (Mat)prob; PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; PetscFunctionBegin; ierr = GetDataID ("lapack","lambda-max-by-real-part-re",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataGetReal ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) rv->r = v; else { ierr = eigenvaluecomp(A); CHKERRQ(ierr); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-real-part-re by now"); } PetscFunctionReturn(0); }
static PetscErrorCode MinEVbyMagIm | ( | AnaModNumericalProblem | prob, |
AnalysisItem * | rv, | ||
int * | lv, | ||
PetscBool * | flg | ||
) | [static] |
Definition at line 297 of file lapack.c.
References eigenvaluecomp(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.
Referenced by RegisterLapackModules().
{ Mat A = (Mat)prob; PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; PetscFunctionBegin; ierr = GetDataID ("lapack","lambda-min-by-magnitude-im",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataGetReal ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) rv->r = v; else { ierr = eigenvaluecomp(A); CHKERRQ(ierr); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-min-by-magnitude-im by now"); } PetscFunctionReturn(0); }
static PetscErrorCode MinEVbyMagRe | ( | AnaModNumericalProblem | prob, |
AnalysisItem * | rv, | ||
int * | lv, | ||
PetscBool * | flg | ||
) | [static] |
Definition at line 272 of file lapack.c.
References eigenvaluecomp(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.
Referenced by RegisterLapackModules().
{ Mat A = (Mat)prob; PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; PetscFunctionBegin; ierr = GetDataID ("lapack","lambda-min-by-magnitude-re",&id,&has); CHKERRQ(ierr); HASTOEXIST(has); ierr = PetscObjectComposedDataGetReal ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) rv->r = v; else { ierr = eigenvaluecomp(A); CHKERRQ(ierr); ierr = PetscObjectComposedDataGetScalar ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); if (*flg) { rv->r = v; } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-min-by-magnitude-re by now"); } PetscFunctionReturn(0); }
PetscErrorCode RegisterLapackModules | ( | ) |
Declare norm-like modules that can be performed with lapack calculations.
Definition at line 424 of file lapack.c.
References ANALYSISDOUBLE, Departure(), MaxEVbyImIm(), MaxEVbyImRe(), MaxEVbyMagIm(), MaxEVbyMagRe(), MaxEVbyRealIm(), MaxEVbyRealRe(), MinEVbyMagIm(), MinEVbyMagRe(), and RegisterModule().
Referenced by AnaModRegisterSalsaModules(), and AnaModRegisterStandardModules().
{ PetscErrorCode ierr; PetscFunctionBegin; ierr = RegisterModule ("lapack","departure",ANALYSISDOUBLE,&Departure); CHKERRQ(ierr); ierr = RegisterModule ("lapack","lambda-max-by-magnitude-re", ANALYSISDOUBLE,&MaxEVbyMagRe); CHKERRQ(ierr); ierr = RegisterModule ("lapack","lambda-max-by-magnitude-im", ANALYSISDOUBLE,&MaxEVbyMagIm); CHKERRQ(ierr); ierr = RegisterModule ("lapack","lambda-min-by-magnitude-re", ANALYSISDOUBLE,&MinEVbyMagRe); CHKERRQ(ierr); ierr = RegisterModule ("lapack","lambda-min-by-magnitude-im", ANALYSISDOUBLE,&MinEVbyMagIm); CHKERRQ(ierr); ierr = RegisterModule ("lapack","lambda-max-by-real-part-re", ANALYSISDOUBLE,&MaxEVbyRealRe); CHKERRQ(ierr); ierr = RegisterModule ("lapack","lambda-max-by-real-part-im", ANALYSISDOUBLE,&MaxEVbyRealIm); CHKERRQ(ierr); ierr = RegisterModule ("lapack","lambda-max-by-im-part-re", ANALYSISDOUBLE,&MaxEVbyImRe); CHKERRQ(ierr); ierr = RegisterModule ("lapack","lambda-max-by-im-part-im", ANALYSISDOUBLE,&MaxEVbyImIm); CHKERRQ(ierr); PetscFunctionReturn(0); }