LCPQ/DEHam
2
0
Fork 0
mirror of https://github.com/LCPQ/DEHam synced 2025-01-03 01:55:55 +01:00
Code Issues Releases Wiki Activity

many new features:

Browse Source 
1. added wf printing support
2. fixed problem with factorial calculation
This commit is contained in:
vijay gopal chilkuri 2018-01-29 16:21:56 +01:00
parent aeabef5ded
commit c8e6e61603
9 changed files with 83 additions and 87 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/elem_diag.irp.f
View File

@ -28,6 +28,10 @@
! if(yw)write(6,*)iaa,'diag,v1' ! if(yw)write(6,*)iaa,'diag,v1'
! endif ! endif
enddo enddo
do i=1, natom
if(deter(i).ne.3) xmatd = xmatd + E(i)
enddo
xmatd = xmatd - E(natom+1)
!-----stockage de l element diag !-----stockage de l element diag

58
src/ex1.c
View File

@ -33,10 +33,7 @@ int main(int argc,char **argv)
PetscInt i,Istart,Iend,col[natomax],maxit,its,nconv,countcol; PetscInt i,Istart,Iend,col[natomax],maxit,its,nconv,countcol;
PetscInt nev, ncv, mpd; PetscInt nev, ncv, mpd;
PetscLogDouble t1,t2,tt1,tt2; PetscLogDouble t1,t2,tt1,tt2;
//PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
PetscErrorCode ierr; PetscErrorCode ierr;
//PetscScalar eigr;
//PetscScalar eigi;
int mpiid; int mpiid;
char const* const fileName = argv[1]; char const* const fileName = argv[1];
@ -44,16 +41,14 @@ int main(int argc,char **argv)
Data getdata; Data getdata;
PetscInt nlocal; PetscInt nlocal;
/* gather the input data */
Data_new(file, &getdata); Data_new(file, &getdata);
getdata.n = get_ntot(getdata.FAM1, getdata.natom, getdata.isz, getdata.ntrou, getdata.fix_trou1, getdata.fix_trou2); getdata.n = get_ntot(getdata.FAM1, getdata.natom, getdata.isz, getdata.ntrou, getdata.fix_trou1, getdata.fix_trou2);
nlocal = getdata.n/getdata.npar; nlocal = getdata.n/getdata.npar;
//printf("n=%ld\t nsites=%d\t nnz=%ld\t npar=%ld\t ntrou=%ld\t isz=%ld\n",getdata.n,getdata.natom, getdata.nnz,getdata.npar,getdata.ntrou,getdata.isz);
PetscScalar *valxr; PetscScalar *valxr;
PetscInt indxr[nlocal]; PetscInt indxr[nlocal];
//Vec Vr,Vi;
char filename[PETSC_MAX_PATH_LEN]="FIL666"; char filename[PETSC_MAX_PATH_LEN]="FIL666";
PetscViewer viewer; PetscViewer viewer;
PetscBool ishermitian; PetscBool ishermitian;
@ -98,30 +93,18 @@ int main(int argc,char **argv)
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,10*getdata.natom,NULL,10*getdata.natom,NULL,&A);CHKERRQ(ierr); ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,10*getdata.natom,NULL,10*getdata.natom,NULL,&A);CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(A,10*getdata.natom,NULL,10*getdata.natom,NULL);CHKERRQ(ierr); ierr = MatMPIAIJSetPreallocation(A,10*getdata.natom,NULL,10*getdata.natom,NULL);CHKERRQ(ierr);
//ierr = MatSetFromOptions(A);CHKERRQ(ierr);
//ierr = MatSetUp(A);CHKERRQ(ierr);
MPI_Comm_rank(MPI_COMM_WORLD,&mpiid); MPI_Comm_rank(MPI_COMM_WORLD,&mpiid);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
ierr = PetscTime(&tt1);CHKERRQ(ierr); ierr = PetscTime(&tt1);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," start: %d end: %d | %d\n",Istart, Iend,getdata.natom);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD," start: %d end: %d \n",Istart, Iend);CHKERRQ(ierr);
// Iend = range[mpiid];
// if(mpiid==0){
// Istart = 0;
// }
// else{
// Istart = range[mpiid-1];
// }
for (i=Istart; i<Iend; i+=getdata.nnz) { for (i=Istart; i<Iend; i+=getdata.nnz) {
tcountcol2=0; tcountcol2=0;
for(kk=0;kk<getdata.nnz;kk++){ for(kk=0;kk<getdata.nnz;kk++){
tcountcol[kk]=0; tcountcol[kk]=0;
} }
iii=i+1; iii=i+1;
// if(i%getdata.npar == 0 && mpiid==0){
// ierr = PetscTime(&t1);CHKERRQ(ierr);
// }
unit_l1_( unit_l1_(
getdata.l1, getdata.l1,
getdata.l2, getdata.l2,
@ -131,6 +114,7 @@ int main(int argc,char **argv)
getdata.xjjxy, getdata.xjjxy,
getdata.xjjz , getdata.xjjz ,
getdata.xtt , getdata.xtt ,
getdata.E ,
tcountcol, tcountcol,
&getdata.ntrou, &getdata.ntrou,
&getdata.isz, &getdata.isz,
@ -147,7 +131,6 @@ int main(int argc,char **argv)
for(kk=0;kk<tcountcol[ll]+1;kk++){ for(kk=0;kk<tcountcol[ll]+1;kk++){
value[kk] = val[kk+tcountcol2]; value[kk] = val[kk+tcountcol2];
col[kk] = tcol[kk+tcountcol2]-1; col[kk] = tcol[kk+tcountcol2]-1;
// PetscPrintf(PETSC_COMM_WORLD,"value = %f col = %d\n",value[kk],col[kk]);
} }
for(kk=tcountcol2+tcountcol[ll]+1;kk<natomax;kk++){ for(kk=tcountcol2+tcountcol[ll]+1;kk<natomax;kk++){
value[kk] = 0.0; value[kk] = 0.0;
@ -169,8 +152,6 @@ int main(int argc,char **argv)
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = PetscTime(&tt2);CHKERRQ(ierr); ierr = PetscTime(&tt2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to assemble the matrix: %f\n",tt2-tt1);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to assemble the matrix: %f\n",tt2-tt1);CHKERRQ(ierr);
//ierr = MatGetVecs(A,NULL,&xr);CHKERRQ(ierr);
//ierr = MatGetVecs(A,NULL,&xi);CHKERRQ(ierr);
ierr = MatCreateVecs(A,NULL,&xr);CHKERRQ(ierr); ierr = MatCreateVecs(A,NULL,&xr);CHKERRQ(ierr);
ierr = MatCreateVecs(A,NULL,&xi);CHKERRQ(ierr); ierr = MatCreateVecs(A,NULL,&xi);CHKERRQ(ierr);
@ -208,9 +189,8 @@ int main(int argc,char **argv)
/* /*
Save eigenvectors, if == ested Save eigenvectors, if == ested
*/ */
//PetscOptionsGetString(NULL,NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs);
EPSGetConverged(eps,&nconv); EPSGetConverged(eps,&nconv);
if (0) { if (getdata.print_wf) {
PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer); PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB); PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU); PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU);
@ -255,10 +235,6 @@ int main(int argc,char **argv)
norm4 = 0.0; norm4 = 0.0;
// for (ii=Istart; ii<Iend; ii+=1) {
// indxr[ii-Istart] = ii;
// }
// ierr = PetscTime(&tt1);CHKERRQ(ierr); // ierr = PetscTime(&tt1);CHKERRQ(ierr);
// ierr = VecGetArray(xr, &valxr);CHKERRQ(ierr); // ierr = VecGetArray(xr, &valxr);CHKERRQ(ierr);
VecScatterCreateToAll(xr,&scatter,&vec2); VecScatterCreateToAll(xr,&scatter,&vec2);
@ -336,34 +312,12 @@ int main(int argc,char **argv)
W3=weight3fin/normfin2; W3=weight3fin/normfin2;
// W3=weight3fin; // W3=weight3fin;
} }
// ierr = PetscTime(&tt2);CHKERRQ(ierr);
// ierr = PetscPrintf(PETSC_COMM_WORLD," Time used to calc par S^2: %f\n",tt2-tt1);CHKERRQ(ierr);
//PetscPrintf(PETSC_COMM_WORLD,"\n norm = %18f xymat = %18f S^2 = %18f \n", norm4, xymat4, norm3);
xymatfin = 0.0; xymatfin = 0.0;
normfin = 0.0; normfin = 0.0;
/* sequential version of analyse */
// ierr = PetscTime(&tt1);CHKERRQ(ierr);
// VecScatterCreateToAll(xr,&scatter,&vec2);
// VecScatterBegin(scatter,xr,vec2,INSERT_VALUES,SCATTER_FORWARD);
// VecScatterEnd(scatter,xr,vec2,INSERT_VALUES,SCATTER_FORWARD);
// VecGetArray(vec2,&values);
// if(mpiid == 0){
// Istart = 0;
// Iend = getdata.n;
// analyse_(values, (Iend-Istart), &Istart, &Iend, &xymat, &norm);
// XS=(1.0/2.0)*(-1.0+sqrt(1.0+(4.0*xymatfin/normfin)));
// printf("\n norm = %18f xymat = %18f S^2 = %18f \n", norm, xymat, XS);
// }
// VecRestoreArray(vec2,&values);
// ierr = PetscTime(&tt2);CHKERRQ(ierr);
// PetscPrintf(PETSC_COMM_WORLD,"seq time = %18f\n",tt2-tt1);
/* /*
Compute the relative error associated to each eigenpair * Compute the relative error associated to each eigenpair
*/ */
ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRQ(ierr); ierr = EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);CHKERRQ(ierr);
@ -385,8 +339,6 @@ int main(int argc,char **argv)
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");CHKERRQ(ierr);
} }
//VecScatterDestroy(&scatter);
//VecDestroy(&vec2);
ierr = EPSDestroy(&eps);CHKERRQ(ierr); ierr = EPSDestroy(&eps);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = VecDestroy(&xr);CHKERRQ(ierr); ierr = VecDestroy(&xr);CHKERRQ(ierr);

8
src/get_ntot.c
View File

@ -15,8 +15,8 @@ int get_ntot(_Bool FAM1, int natom, long int isz, long int ntrou, long int fix_t
natom2 = natom; natom2 = natom;
} }
tnt1 = (int)ceil(exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1))))); tnt1 = (lrint)(exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))));
printf("%10.5f | tnt1=%d\n",exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))),tnt1); //printf("%10.5f | tnt1=%d\n",exp(lgamma((double)(natom2+1)) - (lgamma((double)(natom2-ntrou+1)) + lgamma((double)(ntrou+1)))),tnt1);
int nalpha, nbeta; int nalpha, nbeta;
if((((natom-ntrou) + 2*isz) % 2) == 0){ if((((natom-ntrou) + 2*isz) % 2) == 0){
@ -34,7 +34,7 @@ int get_ntot(_Bool FAM1, int natom, long int isz, long int ntrou, long int fix_t
} }
} }
tnt2 = (int)ceil(exp(lgamma((double)(natom-ntrou+1)) - (lgamma((double)(nalpha+1)) + lgamma((double)(nbeta+1))))); tnt2 = (lrint)(exp(lgamma((double)(natom-ntrou+1)) - (lgamma((double)(nalpha+1)) + lgamma((double)(nbeta+1)))));
printf("natom2=%d fix_trou1=%d fix_trou2=%d nalpha=%d nbeta=%d | | %d %d ntot=%d\n",natom2, fix_trou1, fix_trou2, nalpha, nbeta, tnt1, tnt2, tnt1*tnt2); //printf("natom2=%d fix_trou1=%d fix_trou2=%d nalpha=%d nbeta=%d | | %d %d ntot=%d\n",natom2, fix_trou1, fix_trou2, nalpha, nbeta, tnt1, tnt2, tnt1*tnt2);
return tnt1*tnt2; return tnt1*tnt2;
} }

1
src/get_ntot.h
View File

@ -1,3 +1,4 @@
#include <tgmath.h> #include <tgmath.h>
#include <math.h>
int get_ntot(_Bool Fam1, int natom, long int isz, long int ntrou, long int fix_trou1, long int fix_trou2); int get_ntot(_Bool Fam1, int natom, long int isz, long int ntrou, long int fix_trou1, long int fix_trou2);

1
src/provide_clust.irp.f
View File

@ -4,6 +4,7 @@ BEGIN_PROVIDER[integer,l1, (maxlien)]
&BEGIN_PROVIDER[real*8, xtt ,(maxlien)] &BEGIN_PROVIDER[real*8, xtt ,(maxlien)]
&BEGIN_PROVIDER[real*8, xjjz ,(maxlien)] &BEGIN_PROVIDER[real*8, xjjz ,(maxlien)]
&BEGIN_PROVIDER[real*8, xjjxy,(maxlien)] &BEGIN_PROVIDER[real*8, xjjxy,(maxlien)]
&BEGIN_PROVIDER[real*8, E,(maxlien)]
&BEGIN_PROVIDER[integer, ntrou] &BEGIN_PROVIDER[integer, ntrou]
&BEGIN_PROVIDER[integer, isz] &BEGIN_PROVIDER[integer, isz]
&BEGIN_PROVIDER[logical*1, FAM1] &BEGIN_PROVIDER[logical*1, FAM1]

77
src/read2.c
View File

@ -11,9 +11,9 @@ void Data_new(FILE* file, Data* dat) {
while (fgets(line, sizeof(line), file)) { while (fgets(line, sizeof(line), file)) {
/* note that fgets don't strip the terminating \n, checking its /* note that fgets doesn't strip the terminating \n, checking its
presence would allow to handle lines longer that sizeof(line) */ presence would allow to handle lines longer than sizeof(line) */
if (count != 28){ if (count != 30){
count++; count++;
switch(count){ switch(count){
case 1: case 1:
@ -178,54 +178,84 @@ void Data_new(FILE* file, Data* dat) {
} }
} }
break; break;
case 13: case 13:
dat->nroots=atol(line); arrayIdx=0;
for (token = strtok(line, delim); token != NULL; token = strtok(NULL, delim))
{
double val;
char *unconverted;
/**
* Convert the next token to a float value
*/
val = strtof(token, &unconverted);
if (!isspace(*unconverted) && *unconverted != 0)
{
/**
* Bad input string. Again, we just bail.
*/
fprintf(stderr, "\"%s\" is not a valid floating-point number\n", token);
break;
}
else
{
dat->E[arrayIdx++] = val;
}
}
break; break;
case 14: case 14:
dat->s21a1=atol(line); dat->nroots=atol(line);
break; break;
case 15: case 15:
dat->s21a2=atol(line); dat->s21a1=atol(line);
break; break;
case 16: case 16:
dat->s21b1=atol(line); dat->s21a2=atol(line);
break; break;
case 17: case 17:
dat->s21b2=atol(line); dat->s21b1=atol(line);
break; break;
case 18: case 18:
dat->s22a1=atol(line); dat->s21b2=atol(line);
break; break;
case 19: case 19:
dat->s22a2=atol(line); dat->s22a1=atol(line);
break; break;
case 20: case 20:
dat->s22b1=atol(line); dat->s22a2=atol(line);
break; break;
case 21: case 21:
dat->s22b2=atol(line); dat->s22b1=atol(line);
break; break;
case 22: case 22:
dat->s23a1=atol(line); dat->s22b2=atol(line);
break; break;
case 23: case 23:
dat->s23a2=atol(line); dat->s23a1=atol(line);
break; break;
case 24: case 24:
dat->s23b1=atol(line); dat->s23a2=atol(line);
break; break;
case 25: case 25:
dat->s23b2=atol(line); dat->s23b1=atol(line);
break; break;
case 26: case 26:
dat->postrou=atol(line); dat->s23b2=atol(line);
break; break;
case 27: case 27:
dat->fix_trou1=atol(line); dat->postrou=atol(line);
break; break;
case 28: case 28:
dat->fix_trou1=atol(line);
break;
case 29:
dat->fix_trou2=atol(line); dat->fix_trou2=atol(line);
break; break;
case 30:
dat->print_wf = atol(line);
break;
default:
printf("Done reading file\n");
break;
} /* end of switch */ } /* end of switch */
} /* end of the input file */ } /* end of the input file */
@ -236,10 +266,13 @@ void Data_new(FILE* file, Data* dat) {
} }
_Bool to_bool(const char* str) { _Bool to_bool(const char* str) {
PetscBool strflg; PetscBool strflg=PETSC_FALSE;
PetscStrcmp("true\n",str, &strflg); PetscStrcmp("true\n",str, &strflg);
if(!strflg) PetscStrcmp("True\n",str, &strflg); if(strflg != PETSC_TRUE) PetscStrcmp("True\n",str, &strflg);
if(!strflg) PetscStrcmp("TRUE\n",str, &strflg); if(strflg != PETSC_TRUE) PetscStrcmp("TRUE\n",str, &strflg);
if(strflg != PETSC_TRUE) PetscStrcmp("true",str, &strflg);
if(strflg != PETSC_TRUE) PetscStrcmp("True",str, &strflg);
if(strflg != PETSC_TRUE) PetscStrcmp("TRUE",str, &strflg);
return strflg; return strflg;
} }

16
src/read2.h
View File

@ -12,13 +12,14 @@ typedef struct {
PetscInt n; PetscInt n;
long int nnz,npar; long int nnz,npar;
long int ntrou,isz; long int ntrou,isz;
_Bool FAM1; _Bool FAM1;
long int l1[700]; long int l1[900];
long int l2[700]; long int l2[900];
long int ktyp[700]; long int ktyp[900];
double xjjz[700]; double xjjz[900];
double xjjxy[700]; double xjjxy[900];
double xtt[700]; double xtt[900];
double E[900];
long int nroots; long int nroots;
int natom; int natom;
int s21a1; int s21a1;
@ -36,6 +37,7 @@ typedef struct {
int postrou; int postrou;
long int fix_trou1; long int fix_trou1;
long int fix_trou2; long int fix_trou2;
long int print_wf;
} Data ; } Data ;

1
src/stimsyr.h
View File

@ -9,6 +9,7 @@ void unit_l1_(
double *, double *,
double *, double *,
double *, double *,
double *,
long int *, long int *,
long int *, long int *,
long int *, long int *,

4
src/unit_l1.irp.f
View File

@ -6,6 +6,7 @@
txjjxy, & txjjxy, &
txjjz , & txjjz , &
txtt , & txtt , &
tE , &
tcountcol, & tcountcol, &
tntrou, & tntrou, &
tisz, & tisz, &
@ -24,7 +25,7 @@
integer(C_SIZE_T),INTENT(INOUT)::tcol(maxlien) integer(C_SIZE_T),INTENT(INOUT)::tcol(maxlien)
integer(C_SIZE_T),INTENT(INOUT),dimension(tnrows)::tcountcol integer(C_SIZE_T),INTENT(INOUT),dimension(tnrows)::tcountcol
integer(C_SIZE_T),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien) integer(C_SIZE_T),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien)
real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien) real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien), tE(maxlien)
nrows = tnrows nrows = tnrows
provide nrows provide nrows
@ -35,6 +36,7 @@
xtt(i) = txtt(i) xtt(i) = txtt(i)
xjjxy(i) = txjjxy(i) xjjxy(i) = txjjxy(i)
xjjz (i) = txjjz (i) xjjz (i) = txjjz (i)
E (i) = tE (i)
enddo enddo
ntrou = tntrou ntrou = tntrou
isz = tisz isz = tisz