Skip to content

Commit 6331ad4

Browse files
js1019js1019
authored
Merge pull request #1 from js1019/cHermitian
Update the Hermitian version
2 parents 7bbf93b + 8fb423b commit 6331ad4

18 files changed

+2168
-10
lines changed

‎EXTERNAL/makefile

+2-1
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,9 @@ else
77
#include DIRECTSOL/makefile_mumps.in
88
#OBJ = DIRECTSOL/pevsl_mumps.o ITERSOL/chebiter.o
99
#INCLUDES = -I../INC -I. $(MUMPS_INC)
10+
#include DIRECTSOL/makefile_mumps.in
1011
OBJ = ITERSOL/chebiter.o
11-
INCLUDES = -I../INC -I.
12+
INCLUDES = -I../INC -I.
1213
endif
1314

1415
%.o : %.c

‎FORTRAN/pevsl_f90.c

+233
Original file line numberDiff line numberDiff line change
@@ -81,6 +81,7 @@ void PEVSL_FORT(pevsl_setamv)(uintptr_t *pevslf90, void *func, void *data) {
8181
pEVSL_SetAMatvec(pevsl, (MVFunc) func, data);
8282
}
8383

84+
8485
/** @brief Fortran interface for pEVSL_SetBMatvec
8586
* @param[in] pevslf90 pevsl pointer
8687
* @param[in] func function pointer
@@ -103,6 +104,46 @@ void PEVSL_FORT(pevsl_setbsol)(uintptr_t *pevslf90, void *func, void *data) {
103104
pEVSL_SetBSol(pevsl, (SVFunc) func, data);
104105
}
105106

107+
108+
/**JS 01/20/19 for complex Hermitian Av
109+
* @brief Fortran interface for pEVSL_SetAMatvec
110+
* @param[in] pevslf90 pevsl pointer
111+
* @param[in] func function pointer
112+
* @param[in] data associated data
113+
*/
114+
void PEVSL_FORT(pevsl_setzamv)(uintptr_t *pevslf90, void *func, void *data) {
115+
116+
/* cast pointer */
117+
pevsl_Data *pevsl = (pevsl_Data *) (*pevslf90);
118+
119+
pEVSL_SetZAMatvec(pevsl, (ZMVFunc) func, data);
120+
}
121+
122+
/** @brief Fortran interface for pEVSL_SetBMatvec
123+
* @param[in] pevslf90 pevsl pointer
124+
* @param[in] func function pointer
125+
* @param[in] data associated data
126+
*/
127+
void PEVSL_FORT(pevsl_setzbmv)(uintptr_t *pevslf90, void *func, void *data) {
128+
129+
/* cast pointer */
130+
pevsl_Data *pevsl = (pevsl_Data *) (*pevslf90);
131+
132+
pEVSL_SetZBMatvec(pevsl, (ZMVFunc) func, data);
133+
}
134+
135+
/** @brief Fortran interface for SetBsol */
136+
void PEVSL_FORT(pevsl_setzbsol)(uintptr_t *pevslf90, void *func, void *data) {
137+
138+
/* cast pointer */
139+
pevsl_Data *pevsl = (pevsl_Data *) (*pevslf90);
140+
141+
pEVSL_SetZBSol(pevsl, (ZSVFunc) func, data);
142+
}
143+
144+
145+
146+
106147
/** @brief Fortran interface for SetStdEig */
107148
void PEVSL_FORT(pevsl_set_stdeig)(uintptr_t *pevslf90) {
108149

@@ -182,6 +223,8 @@ void PEVSL_FORT(pevsl_bsv)(uintptr_t *pevslf90, double *b, double *x) {
182223
}
183224

184225

226+
227+
185228
/** @brief Fortran interface for evsl_lanbounds
186229
* @param[in] pevslf90 pEVSL pointer
187230
* @param[in] mlan Krylov dimension
@@ -207,10 +250,46 @@ void PEVSL_FORT(pevsl_lanbounds)(uintptr_t *pevslf90, int *mlan, int *nsteps, do
207250
pEVSL_ParvecRand(&vinit);
208251
/*------------------- Lanczos Bounds */
209252
pEVSL_LanTrbounds(pevsl, *mlan, *nsteps, 1e-8, &vinit, 1, lmin, lmax, NULL);
253+
//pEVSL_LanTrbounds(pevsl, *mlan, *nsteps, 1e-8, &vinit, 1, lmin, lmax, stdout);
210254

211255
pEVSL_ParvecFree(&vinit);
212256
}
213257

258+
259+
/** JS 01/20/19
260+
* @brief Fortran interface for evsl_lanbounds
261+
* @param[in] pevslf90 pEVSL pointer
262+
* @param[in] mlan Krylov dimension
263+
* @param[in] nsteps number of steps
264+
* @param[in] tol stopping tol
265+
* @param[out] lmin lower bound
266+
* @param[out] lmax upper bound
267+
* */
268+
void PEVSL_FORT(pevsl_zlanbounds)(uintptr_t *pevslf90, int *mlan, int *nsteps, double *tol,
269+
double *lmin, double *lmax) {
270+
int N, n, nfirst;
271+
pevsl_Parvec vrinit, viinit;
272+
273+
/* cast pointer */
274+
pevsl_Data *pevsl = (pevsl_Data *) (*pevslf90);
275+
276+
N = pevsl->N;
277+
n = pevsl->n;
278+
nfirst = pevsl->nfirst;
279+
280+
/*------------------- Create parallel vector: random initial guess */
281+
pEVSL_ParvecCreate(N, n, nfirst, pevsl->comm, &vrinit);
282+
pEVSL_ParvecRand_split(&vrinit);
283+
pEVSL_ParvecCreate(N, n, nfirst, pevsl->comm, &viinit);
284+
pEVSL_ParvecRand_split(&viinit);
285+
/*------------------- Lanczos Bounds */
286+
pEVSL_ZLanTrbounds(pevsl, *mlan, *nsteps, 1e-8, &vrinit, &viinit, 1, lmin, lmax, stdout);
287+
288+
pEVSL_ParvecFree(&vrinit);
289+
pEVSL_ParvecFree(&viinit);
290+
}
291+
292+
214293
/** @brief Fortran interface for find_pol
215294
* @param[in] xintv Interval of interest
216295
* @param[in] thresh_int Interior threshold
@@ -291,6 +370,96 @@ void PEVSL_FORT(pevsl_cheblannr)(uintptr_t *pevslf90, double *xintv, int *max_it
291370
pevsl->evec_computed = Y;
292371
}
293372

373+
374+
/** added by JS 020619
375+
* @brief Fortran interface for ChebLanNr
376+
* the results will be saved in the internal variables
377+
*/
378+
void PEVSL_FORT(pevsl_zcheblannr)(uintptr_t *pevslf90, double *xintv, int *max_its, double *tol,
379+
uintptr_t *polf90) {
380+
381+
int N, n, nfirst, nev2, ierr;
382+
double *lam, *res;
383+
pevsl_Parvecs *Yr, *Yi;
384+
FILE *fstats = stdout;
385+
pevsl_Parvec vrinit, viinit;
386+
387+
/* cast pointer */
388+
pevsl_Data *pevsl = (pevsl_Data *) (*pevslf90);
389+
390+
N = pevsl->N;
391+
n = pevsl->n;
392+
nfirst = pevsl->nfirst;
393+
/*-------------------- zero out stats */
394+
pEVSL_StatsReset(pevsl);
395+
/*------------------- Create parallel vector: random initial guess */
396+
pEVSL_ParvecCreate(N, n, nfirst, pevsl->comm, &vrinit);
397+
pEVSL_ParvecCreate(N, n, nfirst, pevsl->comm, &viinit);
398+
pEVSL_ParvecRand(&vrinit);
399+
pEVSL_ParvecRand(&viinit);
400+
/* cast pointer of pol*/
401+
pevsl_Polparams *pol = (pevsl_Polparams *) (*polf90);
402+
/* call ChebLanNr */
403+
ierr = pEVSL_ZChebLanNr(pevsl, xintv, *max_its, *tol, &vrinit, &viinit, pol, &nev2, &lam, &Yr, &Yi, &res, fstats);
404+
405+
if (ierr) {
406+
printf("ZChebLanNr error %d\n", ierr);
407+
}
408+
409+
pEVSL_ParvecFree(&vrinit);
410+
pEVSL_ParvecFree(&viinit);
411+
/*--------------------- print stats */
412+
pEVSL_StatsPrint(pevsl, fstats);
413+
414+
if (res) {
415+
PEVSL_FREE(res);
416+
}
417+
/* save pointers to the global variables */
418+
pevsl->nev_computed = nev2;
419+
pevsl->eval_computed = lam;
420+
pevsl->evec_computed = Yr;
421+
pevsl->evec_imag_computed = Yi;
422+
423+
}
424+
425+
426+
/** JS 051919 for Lanczos vectors
427+
* @brief Fortran interface for ChebLanNr
428+
* the results will be saved in the internal variables
429+
*/
430+
void PEVSL_FORT(pevsl_lanvectors)(uintptr_t *pevslf90, double *xintv, int *max_its, double *tol,
431+
uintptr_t *polf90) {
432+
int N, n, nfirst, nev2, ierr;
433+
pevsl_Parvecs *Y;
434+
FILE *fstats = stdout;
435+
pevsl_Parvec vinit;
436+
437+
/* cast pointer */
438+
pevsl_Data *pevsl = (pevsl_Data *) (*pevslf90);
439+
440+
N = pevsl->N;
441+
n = pevsl->n;
442+
nfirst = pevsl->nfirst;
443+
/*-------------------- zero out stats */
444+
pEVSL_StatsReset(pevsl);
445+
/*------------------- Create parallel vector: random initial guess */
446+
pEVSL_ParvecCreate(N, n, nfirst, pevsl->comm, &vinit);
447+
pEVSL_ParvecRand(&vinit);
448+
/* cast pointer of pol*/
449+
pevsl_Polparams *pol = (pevsl_Polparams *) (*polf90);
450+
/* lanvectors.c */
451+
ierr = pEVSL_Lanvectors(pevsl, xintv, *max_its, *tol, &vinit, pol, &nev2, &Y, fstats);
452+
453+
if (ierr) {
454+
printf("lanvectors error %d\n", ierr);
455+
}
456+
/* save pointers to the global variables */
457+
pevsl->nev_computed = nev2;
458+
pevsl->evec_computed = Y;
459+
//printf("lanvectors done %d\n", ierr);
460+
}
461+
462+
294463
/** @brief Get the number of last computed eigenvalues
295464
*/
296465
void PEVSL_FORT(pevsl_get_nev)(uintptr_t *pevslf90, int *nev) {
@@ -330,6 +499,70 @@ void PEVSL_FORT(pevsl_copy_result)(uintptr_t *pevslf90, double *val, double *vec
330499
PEVSL_FREE(pevsl->evec_computed);
331500
}
332501

502+
503+
/** @brief copy the computed eigenvalues and vectors
504+
* @warning after this call the internal saved results will be freed
505+
* @param[in, out] pevslf90 Data to be cast to pevsl_Data
506+
* @param[in] ld leading dimension of vec [ld >= pevsl.n]
507+
* @param[out] vec Lanczos vectors output
508+
*/
509+
void PEVSL_FORT(pevsl_copy_vectors)(uintptr_t *pevslf90, double *vec, int *ld) {
510+
int i;
511+
/* cast pointer */
512+
pevsl_Data *pevsl = (pevsl_Data *) (*pevslf90);
513+
514+
515+
/* copy eigenvectors */
516+
for (i=0; i<pevsl->nev_computed; i++) {
517+
double *dest = vec + i * (*ld);
518+
double *src = pevsl->evec_computed->data + i * pevsl->evec_computed->ld;
519+
memcpy(dest, src, pevsl->n*sizeof(double));
520+
}
521+
522+
/* reset pointers */
523+
pevsl->nev_computed = 0;
524+
pEVSL_ParvecsFree(pevsl->evec_computed);
525+
PEVSL_FREE(pevsl->evec_computed);
526+
}
527+
528+
529+
530+
/* added JS 020619 for complex Hermitian
531+
* @brief copy the computed eigenvalues and vectors
532+
* @warning after this call the internal saved results will be freed
533+
* @param[in, out] pevslf90 Data to be cast to pevsl_Data
534+
* @param[in] ld leading dimension of vec [ld >= pevsl.n]
535+
* @param[out] val Eigenvalue output
536+
* @param[out] vecr, veci Eigenvector output
537+
*/
538+
void PEVSL_FORT(pevsl_copy_zresult)(uintptr_t *pevslf90, double *val, double *vecr, double *veci, int *ld) {
539+
540+
int i;
541+
/* cast pointer */
542+
pevsl_Data *pevsl = (pevsl_Data *) (*pevslf90);
543+
/* copy eigenvalues */
544+
memcpy(val, pevsl->eval_computed, pevsl->nev_computed*sizeof(double));
545+
546+
/* copy eigenvectors */
547+
for (i=0; i<pevsl->nev_computed; i++) {
548+
double *destr = vecr + i * (*ld);
549+
double *srcr = pevsl->evec_computed->data + i * pevsl->evec_computed->ld;
550+
memcpy(destr, srcr, pevsl->n*sizeof(double));
551+
double *desti = veci + i * (*ld);
552+
double *srci = pevsl->evec_imag_computed->data + i*pevsl->evec_imag_computed->ld;
553+
memcpy(desti, srci, pevsl->n*sizeof(double));
554+
}
555+
556+
/* reset pointers */
557+
pevsl->nev_computed = 0;
558+
PEVSL_FREE(pevsl->eval_computed);
559+
pEVSL_ParvecsFree(pevsl->evec_computed);
560+
PEVSL_FREE(pevsl->evec_computed);
561+
562+
pEVSL_ParvecsFree(pevsl->evec_imag_computed);
563+
PEVSL_FREE(pevsl->evec_imag_computed);
564+
}
565+
333566
void PEVSL_FORT(pevsl_setup_chebiter)(double *lmin, double *lmax, int *deg,
334567
uintptr_t *parcsrf90, uintptr_t *chebf90) {
335568
/* cast pointer of the matrix */

‎INC/pevsl.h

+29-1
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
/*
2-
This file contains function prototypes and constant definitions for EVSL
2+
This file contains function prototypes and constant definitions for pEVSL
33
*/
44

55
#ifndef PEVSL_H
@@ -28,6 +28,10 @@ void pEVSL_SpslicerLan(double* xi, double* yi, int n_int, int npts, double* sli)
2828
int pEVSL_LanTrbounds(pevsl_Data *pevsl, int lanm, int maxit, double tol,
2929
pevsl_Parvec *vinit, int bndtype, double *lammin,
3030
double *lammax, FILE *fstats);
31+
/* add JS 010219 complex lantrbnd */
32+
int pEVSL_ZLanTrbounds(pevsl_Data *pevsl, int lanm, int maxit, double tol,
33+
pevsl_Parvec *vrinit, pevsl_Parvec *viinit, int bndtype,
34+
double *lammin, double *lammax, FILE *fstats);
3135
/* lspolapprox.c */
3236
void pEVSL_SetupLSPolSqrt(int max_deg, double tol, double lmin, double lmax,
3337
pevsl_Parcsr *B, void **vdata);
@@ -49,6 +53,10 @@ int pEVSL_SetLTSol (pevsl_Data *pevsl_data, SVFunc func, void *data);
4953
int pEVSL_SetStdEig (pevsl_Data *pevsl_data);
5054
int pEVSL_SetGenEig (pevsl_Data *pevsl_data);
5155
int pEVSL_SetSigmaMult(pevsl_Data *pevsl_data, double mult);
56+
/* add JS 01/20/19 complex mat-vec*/
57+
int pEVSL_SetZAMatvec (pevsl_Data *pevsl_data, ZMVFunc func, void *data);
58+
int pEVSL_SetZBMatvec (pevsl_Data *pevsl_data, ZMVFunc func, void *data);
59+
int pEVSL_SetZBSol (pevsl_Data *pevsl_data, ZSVFunc func, void *data);
5260

5361
/* parcsr.c */
5462
int pEVSL_ParcsrCreate(int nrow, int ncol, int *row_starts, int *col_starts,
@@ -75,6 +83,13 @@ void pEVSL_ParvecNrm2(pevsl_Parvec *x, double *t);
7583
void pEVSL_ParvecCopy(pevsl_Parvec *x, pevsl_Parvec *y);
7684
void pEVSL_ParvecSum(pevsl_Parvec *x, double *t);
7785
void pEVSL_ParvecScal(pevsl_Parvec *x, double t);
86+
/*add JS 01/02/19 */
87+
void pEVSL_ParvecZNrm2(pevsl_Parvec *xr, pevsl_Parvec *xi, double *t);
88+
void pEVSL_ParvecZDot(pevsl_Parvec *xr, pevsl_Parvec *xi, pevsl_Parvec *yr,
89+
pevsl_Parvec *yi, double *tr, double *ti);
90+
void pEVSL_ParvecRand_split(pevsl_Parvec *x);
91+
92+
7893
/* void pEVSL_ParvecAddScalar(pevsl_Parvec *x, double t); */
7994
void pEVSL_ParvecSetScalar(pevsl_Parvec *x, double t);
8095
void pEVSL_ParvecSetZero(pevsl_Parvec *x);
@@ -141,5 +156,18 @@ int pEVSL_ChebLanNr(pevsl_Data *pevsl, double *intv, int maxit, double tol, pevs
141156
pevsl_Polparams *pol, int *nevOut, double **lamo, pevsl_Parvecs **Wo,
142157
double **reso, FILE *fstats);
143158

159+
/* add JS 020619 for complex Hermitian cases*/
160+
int pEVSL_ZChebLanNr(pevsl_Data *pevsl, double *intv, int maxit, double tol,
161+
pevsl_Parvec *vrinit, pevsl_Parvec *viinit, pevsl_Polparams *pol, int *nevOut,
162+
double **lamo, pevsl_Parvecs **Wor, pevsl_Parvecs **Woi,
163+
double **reso, FILE *fstats);
164+
165+
166+
/* add JS 051919 for lanczos vectors */
167+
/* lanvectors.c */
168+
int pEVSL_Lanvectors(pevsl_Data *pevsl, double *intv, int maxit, double tol,
169+
pevsl_Parvec *vinit, pevsl_Polparams *pol, int *nevOut,
170+
pevsl_Parvecs **Wo, FILE *fstats);
171+
144172
#endif
145173

‎INC/pevsl_blaslapack.h

+3
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,9 @@ void DHSEQR(char* jobz,char* compz,int* n,int* ilo,int* ihi,double* h,int* ldh,d
4949
double* z,int* ldz,double* work, int* lwork,int* info);
5050
void ZGESV(int *n, int *nrow, complex double * A, int* m, int* ipiv, complex double *rhs, int* k, int* INFO);
5151

52+
/* add complex Hermitian eigensolver */
53+
54+
5255
#endif
5356

5457
/* USE SQRT(DDOT()) INSTEAD OF DNRM2, FOR BETTER PERFORMANCE */

0 commit comments

Comments
 (0)