From 79bdcfa221a96adaea0bc7daf56507c9f5f5e465 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Mon, 19 Dec 2016 23:26:16 +0100 Subject: [PATCH] initial commit --- adr.irp.f | 39 +++++ adrfull.irp.f | 54 +++++++ combin.irp.f | 28 ++++ conv.irp.f | 20 +++ desort.irp.f | 32 ++++ deter.irp.f | 9 ++ elem_diag.irp.f | 43 ++++++ ex1.c | 212 +++++++++++++++++++++++++ extra_diag.irp.f | 135 ++++++++++++++++ foundet.irp.f | 16 ++ gamma.irp.f | 100 ++++++++++++ getdet.irp.f | 69 +++++++++ main.irp.f | 83 ++++++++++ natom.irp.f | 367 ++++++++++++++++++++++++++++++++++++++++++++ natomax.irp.f | 22 +++ nnk.irp.f | 9 ++ nt1.irp.f | 13 ++ nt2.irp.f | 10 ++ prodcol.irp.f | 8 + prodval.irp.f | 4 + provide_clust.irp.f | 26 ++++ pstartend.irp.f | 6 + rank.irp.f | 15 ++ read2.c | 211 +++++++++++++++++++++++++ read2.h | 20 +++ searchdet.irp.f | 86 +++++++++++ searchdetfull.irp.f | 140 +++++++++++++++++ sort.irp.f | 32 ++++ stimsyr.h | 15 ++ unit_FIL44.irp.f | 93 +++++++++++ unit_l1.irp.f | 45 ++++++ ylogic.irp.f | 94 ++++++++++++ 32 files changed, 2056 insertions(+) create mode 100644 adr.irp.f create mode 100644 adrfull.irp.f create mode 100644 combin.irp.f create mode 100644 conv.irp.f create mode 100644 desort.irp.f create mode 100644 deter.irp.f create mode 100644 elem_diag.irp.f create mode 100644 ex1.c create mode 100644 extra_diag.irp.f create mode 100644 foundet.irp.f create mode 100644 gamma.irp.f create mode 100644 getdet.irp.f create mode 100644 main.irp.f create mode 100644 natom.irp.f create mode 100644 natomax.irp.f create mode 100644 nnk.irp.f create mode 100644 nt1.irp.f create mode 100644 nt2.irp.f create mode 100644 prodcol.irp.f create mode 100644 prodval.irp.f create mode 100644 provide_clust.irp.f create mode 100644 pstartend.irp.f create mode 100644 rank.irp.f create mode 100644 read2.c create mode 100644 read2.h create mode 100644 searchdet.irp.f create mode 100644 searchdetfull.irp.f create mode 100644 sort.irp.f create mode 100644 stimsyr.h create mode 100644 unit_FIL44.irp.f create mode 100644 unit_l1.irp.f create mode 100644 ylogic.irp.f diff --git a/adr.irp.f b/adr.irp.f new file mode 100644 index 0000000..76db98d --- /dev/null +++ b/adr.irp.f @@ -0,0 +1,39 @@ +subroutine adr(ideter,add) + implicit none + BEGIN_DOC + ! this subroutine provides the address of a detrminant + ! given in old format. + ! It searches in a list of generated determinants and + ! matches the given determinant. + END_DOC + integer,INTENT(INOUT)::ideter(natomax) + integer(kind=selected_int_kind(16)),INTENT(INOUT)::add + integer(kind=selected_int_kind(16))::det,deth,addh,detnew + integer::count,i,j + + det=0 + detnew=0 + deth=0 + count=0 + call conv(ideter,det,deth) + Do i=0,natom-1 + if(BTEST(deth,i))then + count=count+1 + endif + if(BTEST(det,i))then + detnew=IBSET(detnew,i-count) + endif + enddo + det=detnew + call searchdet(det,add,deth,addh) + add = add + (nt1-addh)*(nt2) + + +10 FORMAT(B64,I8,F8.2) +15 FORMAT(B64,I8,I8,I8) +11 FORMAT(B64,I3,B64) +12 FORMAT(I5,$) +13 FORMAT(B64,B64) +14 FORMAT(B64,I14) +16 FORMAT(B64,I14,I14) +end diff --git a/adrfull.irp.f b/adrfull.irp.f new file mode 100644 index 0000000..9cd7f65 --- /dev/null +++ b/adrfull.irp.f @@ -0,0 +1,54 @@ +subroutine adrfull() + implicit none + BEGIN_DOC + ! this subroutine provides the address of a detrminant + ! given in old format. + ! It searches in a list of generated determinants and + ! matches the given determinant. + END_DOC + integer,dimension(natomax)::ideter + integer(kind=selected_int_kind(16))::add + integer(kind=selected_int_kind(16))::det,deth,addh,detnew + integer::count,i,j + + det=0 + detnew=0 + deth=0 + do j=1,detfound + detnew=0 + count=0 + ideter=foundet(:,j) + call conv(ideter,det,deth) + Do i=0,natom-1 + if(BTEST(deth,i))then + count=count+1 + endif + if(BTEST(det,i))then + detnew=IBSET(detnew,i-count) + endif + enddo + det=detnew + foundadd(j,1)=det + foundadd(j,3)=j + foundaddh(j,1)=deth + foundaddh(j,3)=j + enddo + call sort() + call searchdetfull() + call desort() + do i=1,detfound + add = foundadd(i,2) + addh = foundaddh(i,2) + add = add + (nt1-addh)*(nt2) + foundetadr(i)=add + enddo + + +10 FORMAT(B64,I8,F8.2) +15 FORMAT(B64,I8,I8,I8) +11 FORMAT(B64,I3,B64) +12 FORMAT(I5,$) +13 FORMAT(B64,B64) +14 FORMAT(B64,I14) +16 FORMAT(B64,I14,I14) +end diff --git a/combin.irp.f b/combin.irp.f new file mode 100644 index 0000000..b96c173 --- /dev/null +++ b/combin.irp.f @@ -0,0 +1,28 @@ + SUBROUTINE COMBIN(NS1,N,M,NT,MD,ND) + INTEGER NS1(MD,ND) + INTEGER :: i,j,MMJ1,J1 + I=1 + DO J=1,M + NS1(j,i)=J + enddo + J1=0 + do while (J1.LT.M) + J1=0 + DO J=1,M + IF(NS1(j,i).EQ.N-M+J) J1=J1+1 + enddo + IF(J1.EQ.M) EXIT + I=I+1 + DO J=1,M + NS1(j,i)=NS1(j,I-1) + enddo + NS1(m-j1,I)=NS1(m-j1,I-1)+1 + IF(J1.EQ.0) CYCLE + MMJ1=M-J1+1 + DO J=MMJ1,M + NS1(j,I)=NS1(j-1,I)+1 + enddo + enddo + NT=I + RETURN + END diff --git a/conv.irp.f b/conv.irp.f new file mode 100644 index 0000000..f296c60 --- /dev/null +++ b/conv.irp.f @@ -0,0 +1,20 @@ +subroutine conv(ideter,det,deth) + implicit none + BEGIN_DOC + ! this routine converts a detrminant in the old + ! format into the new one and returns the determinant. + END_DOC + integer,INTENT(INOUT)::ideter(natomax) + integer(kind=selected_int_kind(16)),INTENT(INOUT)::det + integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth + integer::i + det=0 + deth=0 + do i=1,natom + if(ideter(natom-i+1).eq.2 .and. ideter(natom-i+1).ne.3)then + det=IBSET(det,i-1) + elseif(ideter(natom-i+1).eq.3)then + deth=IBSET(deth,i-1) + endif + enddo +end diff --git a/desort.irp.f b/desort.irp.f new file mode 100644 index 0000000..c4d09b1 --- /dev/null +++ b/desort.irp.f @@ -0,0 +1,32 @@ +subroutine desort() + implicit none + integer::i,j,ord,ordh + integer(kind=selected_int_kind(16))::add,addh,det,deth,addt + + do i=1,detfound-1 + do j=i+1,detfound + if(foundaddh(i,3).gt.foundaddh(j,3))then + deth = foundaddh(i,1) + foundaddh(i,1) = foundaddh(j,1) + foundaddh(j,1) = deth + addh = foundaddh(i,2) + foundaddh(i,2) = foundaddh(j,2) + foundaddh(j,2) = addh + ordh = foundaddh(i,3) + foundaddh(i,3) = foundaddh(j,3) + foundaddh(j,3) = ordh + endif + if(foundadd(i,3).gt.foundadd(j,3))then + det = foundadd(i,1) + foundadd(i,1) = foundadd(j,1) + foundadd(j,1) = det + add = foundadd(i,2) + foundadd(i,2) = foundadd(j,2) + foundadd(j,2) = add + ord = foundadd(i,3) + foundadd(i,3) = foundadd(j,3) + foundadd(j,3) = ord + endif + enddo + enddo +end diff --git a/deter.irp.f b/deter.irp.f new file mode 100644 index 0000000..ade96ba --- /dev/null +++ b/deter.irp.f @@ -0,0 +1,9 @@ +BEGIN_PROVIDER [integer, deter, (natomax)] + + implicit none + BEGIN_DOC + ! provides ideter and iaa + END_DOC + integer ::ideter(natomax) + deter=0 +END_PROVIDER diff --git a/elem_diag.irp.f b/elem_diag.irp.f new file mode 100644 index 0000000..679b0d2 --- /dev/null +++ b/elem_diag.irp.f @@ -0,0 +1,43 @@ + subroutine elem_diag(xmatd) + + implicit none + + integer :: i + real*8 :: xmatd + logical :: yw + +! write(6,*)'in elem_diag' +!--- + +! if(yw)write(6,*)'ideter',(det(iaa,i),i=1,natom) + + yw=.FALSE. + xmatd=0.d0 +! testing code +! xmatd=1000.d0 +! if(.not.redo)write(6,*)'vijayyves' + do i=1,nlientot + if(yalt(i))then +! xmatd=-(xj1+xeneJ(i)*xbJ)+xmatd + xmatd= -(xjz(i))+xmatd + if(yw)write(6,*)xmatd,'xmatd' + if(yw)write(6,*)'xjz',xjz(i) + endif +! if(yrep1(i))then +! xmat=xv1+xmat +! if(yw)write(6,*)iaa,'diag,v1' +! endif + enddo + +!-----stockage de l element diag + +! imat=iaa +! imat3(isto3+1)=imat +! jmat3(isto3+1)=imat +! dmat3(isto3+1)=xmat/2 + + return + end + + + diff --git a/ex1.c b/ex1.c new file mode 100644 index 0000000..b29a2b5 --- /dev/null +++ b/ex1.c @@ -0,0 +1,212 @@ +#include +#include +#include "stimsyr.h" +#include "read2.h" + +#undef __FUNCT__ +#define __FUNCT__ "main" + +int main(int argc,char **argv) +{ + Mat A; /* problem matrix */ + EPS eps; /* eigenproblem solver context */ + EPSType type; + PetscReal error,tol,re,im; + PetscScalar kr,ki,value[700]; + Vec xr,xi; + PetscInt i,Istart,Iend,col[700],maxit,its,nconv,countcol; + PetscInt nev, ncv, mpd; + PetscLogDouble t1,t2,tt1,tt2; +//PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE; + PetscErrorCode ierr; +//PetscScalar eigr; +//PetscScalar eigi; + int mpiid; + + char const* const fileName = argv[1]; + FILE* file = fopen(fileName, "r"); + Data getdata; + + Data_new(file, &getdata); +//printf("n=%ld\t nnz=%ld\t npar=%ld\t ntrou=%ld\t isz=%ld\n",getdata.n,getdata.nnz,getdata.npar,getdata.ntrou,getdata.isz); + + +//Vec Vr,Vi; + char filename[PETSC_MAX_PATH_LEN]="FIL666"; + PetscViewer viewer; + PetscBool ishermitian; + PetscInt kk,ll,iii2; + long int iii; + long int tcountcol2,tcol[700],tcountcol[getdata.nnz]; + double val[700]; + + SlepcInitialize(&argc,&argv,(char*)0,NULL); + ierr = PetscPrintf(PETSC_COMM_WORLD,"\n1-D t-J Eigenproblem, n=%D\n\n",getdata.n);CHKERRQ(ierr); + ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); + ierr = MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,getdata.n,getdata.n,getdata.nnz*getdata.npar,NULL,getdata.nnz*getdata.npar,NULL,&A);CHKERRQ(ierr); + ierr = MatMPIAIJSetPreallocation(A,getdata.nnz*getdata.npar,NULL,getdata.nnz*getdata.npar,NULL);CHKERRQ(ierr); +//ierr = MatSetFromOptions(A);CHKERRQ(ierr); +//ierr = MatSetUp(A);CHKERRQ(ierr); + + MPI_Comm_rank(MPI_COMM_WORLD,&mpiid); + ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); + ierr = PetscTime(&tt1);CHKERRQ(ierr); + + for (i=Istart; i0) { + /* + Display eigenvalues and relative errors + */ + ierr = PetscPrintf(PETSC_COMM_WORLD, + " k ||Ax-kx||/||kx||\n" + " ----------------- ------------------\n");CHKERRQ(ierr); + + for (i=0;i0) { + PetscViewerASCIIOpen(PETSC_COMM_WORLD,filename,&viewer); + PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_MATLAB); + PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_SYMMODU); + EPSIsHermitian(eps,&ishermitian); + for (i=0;i +#include +#include +#include +#include "read2.h" + +void Data_new(FILE* file, Data* dat) { + +//ata* dat = (Data*)malloc(sizeof(Data)); + char line[256]; + char *token; + const char *delim=","; + size_t arrayIdx=0; + int count=0; + + while (fgets(line, sizeof(line), file)) { + + /* note that fgets don't strip the terminating \n, checking its + presence would allow to handle lines longer that sizeof(line) */ + if (count != 11){ + count++; + switch(count){ + case 1: + dat->n=atol(line); + break; + case 2: + dat->nnz=atol(line); + break; + case 3: + dat->npar=atol(line); + break; + case 4: + dat->ntrou=atol(line); + break; + case 5: + dat->isz=atol(line); + break; + case 6: + 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 = strtol(token, &unconverted, 10); + 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->l1[arrayIdx++] = val; + } + } + break; + case 7: + 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 = strtol(token, &unconverted, 10); + 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->l2[arrayIdx++] = val; + } + } + break; + case 8: + 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 = strtol(token, &unconverted, 10); + 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->ktyp[arrayIdx++] = val; + } + } + break; + case 9: + 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->xjjz[arrayIdx++] = val; + } + } + break; + case 10: + 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->xjjxy[arrayIdx++] = val; + } + } + break; + case 11: + 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->xtt[arrayIdx++] = val; + } + } + break; + } /* end of switch */ + + } /* end of the input file */ + + } /* end of loop */ + +//return dat; +} + +/* +int main(int argc, char* argv[]) +{ + char const* const fileName = argv[1]; + FILE* file = fopen(fileName, "r"); + + Data *getdata=Data_new(file); + + fclose(file); + + printf("n=%d,nnz=%d,npar=%d\n",getdata->n,getdata->nnz,getdata->npar); + for(int i=0;i<7;i++){ + printf("l1(%d)=%d l2(%d)=%d l3(%d)=%d\n",i,getdata->l1[i],i,getdata->l2[i],i,getdata->l3[i]); + } + for(int i=0;i<7;i++){ + printf("xjjz(%d)=%f xjjxy(%d)=%f xtt(%d)=%f\n",i,getdata->xjjz[i],i,getdata->xjjxy[i],i,getdata->xtt[i]); + } + + return 0; +} +*/ diff --git a/read2.h b/read2.h new file mode 100644 index 0000000..6c98ab4 --- /dev/null +++ b/read2.h @@ -0,0 +1,20 @@ +#include +#include +#include +#include +#include + +typedef struct { + PetscInt n; + long int nnz,npar; + long int ntrou,isz; + long int l1[700]; + long int l2[700]; + long int ktyp[700]; + double xjjz[700]; + double xjjxy[700]; + double xtt[700]; + +} Data ; + +void Data_new(FILE* , Data* ); diff --git a/searchdet.irp.f b/searchdet.irp.f new file mode 100644 index 0000000..d900984 --- /dev/null +++ b/searchdet.irp.f @@ -0,0 +1,86 @@ +subroutine searchdet(det,add,deth,addh) + BEGIN_DOC + ! this subroutine is at the heart of the idea + ! it will generate all the determinants in a fixed order + ! then find the posistion of the determinant given and + ! return it's position in add. + END_DOC + integer(kind=selected_int_kind(16)),INTENT(INOUT)::det + integer(kind=selected_int_kind(16)),INTENT(INOUT)::add + integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth + integer(kind=selected_int_kind(16)),INTENT(INOUT)::addh + integer(kind=selected_int_kind(16))::dethsh + integer(kind=selected_int_kind(16))::a + integer(kind=selected_int_kind(16))::i + integer::const + i=1 + a=0 + add=0 + const=0 + + If(ntrou.ge.1)then + + const=0 + dethsh = ISHFT(deth,-natom/2) + addh=0 +! i=nt1 + do while (i.le.(nt1)) + if(a.eq.dethsh)then + addh=i-1 + EXIT + endif + + i+=1 + a+=1 + do while(popcnt(a).ne.ntrou) + a+=1 + enddo + enddo + if(a.eq.dethsh .and. addh.eq.0)then + addh=nt1 + endif + + endif + + !C if det=0 then exit + a=0 + i=0 + count=0 + if(a.eq.det)then + add=1 + Return + endif + + do while (i.le.(nt2)) + if(a.eq.det)then + if(a.eq.1)then + add=i + count=-1 + EXIT + else + add=i + count=-1 + EXIT + endif + endif + + i+=1 + a+=1 +!C write(6,16)a,a,i-2 + do while(popcnt(a).ne.nbeta) + a+=1 + enddo + enddo + if(a.eq.det .and. count.ne.-1)then + add=i-1 + endif + + +10 FORMAT(B64,I8,F8.2) +15 FORMAT(B64,I8,I8,I8) +11 FORMAT(B64,I3,B64) +12 FORMAT(I5,$) +13 FORMAT(B64,B64) +14 FORMAT(B64,I8) +16 FORMAT(B64,I8,I8) +end diff --git a/searchdetfull.irp.f b/searchdetfull.irp.f new file mode 100644 index 0000000..362c527 --- /dev/null +++ b/searchdetfull.irp.f @@ -0,0 +1,140 @@ +subroutine searchdetfull() + BEGIN_DOC + ! this subroutine is at the heart of the idea + ! it will generate all the determinants in a fixed order + ! then find the posistion of the determinant given and + ! return it's position in add. + END_DOC +! integer(kind=selected_int_kind(16)),INTENT(INOUT)::foundetadr(maxlien,4) + integer(kind=selected_int_kind(16))::add +! integer(kind=selected_int_kind(16)),INTENT(INOUT)::deth + integer(kind=selected_int_kind(16))::dethsh + integer(kind=selected_int_kind(16))::addh + integer(kind=selected_int_kind(16))::a + integer(kind=selected_int_kind(16))::i + integer::const,count + i=1 + a=0 + add=0 + const=0 + count=1 + If(ntrou.ge.1)then + + const=0 + a=0 + addh=0 + i=1 + dethsh = ISHFT(foundaddh(count,1),-natom/2) + do while (i.le.(2*nt1)) + if(a.eq.dethsh)then + addh=i-1 + foundaddh(count,2)=addh + count+=1 + dethsh = ISHFT(foundaddh(count,1),-natom/2) + do while(count .le. detfound .and. a .eq. dethsh) + foundaddh(count,2)=addh + count+=1 + dethsh = ISHFT(foundaddh(count,1),-natom/2) + enddo + if(count.gt.detfound)then + EXIT + endif + endif + + i+=1 +! const=1 + a+=1 +! do while(popcnt(a).ne.ntrou .or. const==1) + do while(popcnt(a).ne.ntrou) + a+=1 +! const=0 + enddo + enddo + + if(a.eq.foundaddh(count,1))then + addh=i-1 + foundaddh(count,2)=addh + count+=1 + dethsh = ISHFT(foundaddh(count,1),-natom/2) + do while(count .le. detfound .and. a .eq. foundaddh(count,1)) + foundaddh(count,2)=addh + count+=1 + dethsh = ISHFT(foundaddh(count,1),-natom/2) + enddo + endif + + endif + + !C if det=0 then exit + a=0 + i=1 + count=1 + const=0 + if(a.eq.foundadd(count,1))then + add=1 + foundadd(count,2)=add + count+=1 + do while(foundadd(count,1).eq.a .and. count.le.detfound) + foundadd(count,2)=add + count+=1 + enddo + endif + + do while (i.le.(nt2) .and. count .le.detfound) + if(a.eq.foundadd(count,1))then + if(a.eq.1)then + add=i-1 + foundadd(count,2)=add + count+=1 + do while(foundadd(count,1).eq.a .and. count.le.detfound) + foundadd(count,2)=add + count+=1 + enddo + if(count.eq.detfound)then + const=-1 + EXIT + endif + else + add=i-1 + foundadd(count,2)=add + count+=1 + do while(foundadd(count,1).eq.a .and. count.le.detfound) + foundadd(count,2)=add + count+=1 + enddo + if(count.gt.detfound)then + const=-1 + EXIT + endif + endif + endif + + i+=1 + a+=1 +!C write(6,16)a,a,i-2 + do while(popcnt(a).ne.nbeta) + a+=1 + enddo + enddo +! if(a.eq.foundadd(count,1) .and. foundadd(count,1).eq.0)then + if(a.eq.foundadd(count,1) .and. foundadd(count,2) .eq. 0)then + foundadd(count,2) = nt2 + count+=1 + do while(foundadd(count,1).eq.a .and. count.le.detfound) + foundadd(count,2)=nt2 + count+=1 + enddo + endif +! do i=1,detfound +! write(6,16)foundadd(i,1) ,foundadd(i,2),foundadd(i,3) +! write(6,16)foundaddh(i,1),foundaddh(i,2),foundaddh(i,3) +! enddo + +10 FORMAT(B64,I8,F8.2) +15 FORMAT(B64,I8,I8,I8) +11 FORMAT(B64,I3,B64) +12 FORMAT(I5,$) +13 FORMAT(B64,B64) +14 FORMAT(B64,I8) +16 FORMAT(B64,I8,I8) +end diff --git a/sort.irp.f b/sort.irp.f new file mode 100644 index 0000000..1c8c130 --- /dev/null +++ b/sort.irp.f @@ -0,0 +1,32 @@ +subroutine sort() + implicit none + integer::i,j,ord,ordh + integer(kind=selected_int_kind(16))::add,addh,det,deth,addt + + do i=1,detfound-1 + do j=i+1,detfound + if(foundaddh(i,1).gt.foundaddh(j,1))then + deth = foundaddh(i,1) + foundaddh(i,1) = foundaddh(j,1) + foundaddh(j,1) = deth + addh = foundaddh(i,2) + foundaddh(i,2) = foundaddh(j,2) + foundaddh(j,2) = addh + ordh = foundaddh(i,3) + foundaddh(i,3) = foundaddh(j,3) + foundaddh(j,3) = ordh + endif + if(foundadd(i,1).gt.foundadd(j,1))then + det = foundadd(i,1) + foundadd(i,1) = foundadd(j,1) + foundadd(j,1) = det + add = foundadd(i,2) + foundadd(i,2) = foundadd(j,2) + foundadd(j,2) = add + ord = foundadd(i,3) + foundadd(i,3) = foundadd(j,3) + foundadd(j,3) = ord + endif + enddo + enddo +end diff --git a/stimsyr.h b/stimsyr.h new file mode 100644 index 0000000..84dc62f --- /dev/null +++ b/stimsyr.h @@ -0,0 +1,15 @@ +void unit_l1_( + long int *, + long int *, + long int *, + long int *, + long int *, + double *, + double *, + double *, + long int *, + long int *, + long int *, + long int *, + double *); + diff --git a/unit_FIL44.irp.f b/unit_FIL44.irp.f new file mode 100644 index 0000000..647554b --- /dev/null +++ b/unit_FIL44.irp.f @@ -0,0 +1,93 @@ + subroutine unit(tistart,tcountcol,tcol,tval) + + BEGIN_DOC + ! file units for writing + END_DOC + + implicit none + integer :: i,j,k,ia1,ia2,l,m,chcind,chcval,ii,tistart2 + integer :: count,unit_44,unit_33 + integer :: iat,nbtots + integer(kind=selected_int_kind(16))::iaa + integer :: kkio,kkiok,n,nz,cdiag,cexdiag + integer,allocatable ::ideter1(:),ideter2(:),deti(:),detj(:) + integer(kind=selected_int_kind(16)),dimension(maxlien) ::tl1,tl2,tktyp + integer(kind=selected_int_kind(16)),dimension(nrows)::tcountcol + integer(kind=selected_int_kind(16))::tistart + real*8,dimension(maxlien)::tval + integer(kind=selected_int_kind(16)),dimension(maxlien)::tcol + real*8 :: xmat + integer :: ik,imat4,iaa2,iik + integer :: ik1,ik2,jmat4,IC,ikmax,ikmin + real*8 :: dmat4 + logical :: yw + ! BEGIN_DOC + ! provides unit of FIL33 & FIL44 + ! END_DOC + + allocate (ideter2(natomax)) +! allocate (tcol(natomax)) +! allocate (tval(natomax)) + + countcol=0 + unit_44=44 + unit_33=33 + nnk=0 + xmat=0d0 + count=0 + cdiag=1 + cexdiag=0 + tistart2=tistart + + + do i=1,natomax + col(i)=0 + val(i)=0d0 + tval(i)=0d0 + tcol(i)=0 + enddo + tcountcol=0 + countcol=0 + xmat=0d0 + count=0 + +! tistart=tistart + i=1+tistart/nt2 + k=1+mod(tistart , nt2) + +! call getdet(tistart,ideter2) +! deter=ideter2 +! Touch deter +! call adr(deter,iaa) +! call elem_diag(xmat) +! countcol+=1 +! col(countcol)=iaa +! val(countcol)=xmat*1.0d0 + + call extra_diag(tistart) + + tcountcol=countcolfull + do i=1,maxlien + if(col(i).ne.0)then + if(val(i) .ne. 0 .or. col(i).eq.tistart2)then + tcol(i)=col(i) + tval(i)=val(i) + endif + if(col(i).eq.tistart2)then + cexdiag+=1 + elseif(cexdiag .eq. countcolfull(cdiag))then + cexdiag=0 + if(cdiag.lt.nrows)then + cdiag+=1 + endif + tistart2+=1 + else + cexdiag+=1 + endif + endif + enddo +! print *,"tistart =", tistart,"countcol =", countcol,"\n",(tcountcol(i),i=1,nrows) +! print *,"" +! print *,(tcol(i),i=1,maxlien) + deallocate(ideter2) + end diff --git a/unit_l1.irp.f b/unit_l1.irp.f new file mode 100644 index 0000000..5239a49 --- /dev/null +++ b/unit_l1.irp.f @@ -0,0 +1,45 @@ + subroutine unit_l1(tl1,& + tl2, & + tktyp, & + tistart, & + tnrows, & + txjjxy, & + txjjz , & + txtt , & + tcountcol, & + tntrou, & + tisz, & + tcol,tval) + implicit none + integer,INTENT(INOUT)::tistart, tnrows, tntrou, tisz + integer::i + real*8,INTENT(INOUT)::tval(maxlien) + integer(kind=selected_int_kind(16)),INTENT(INOUT)::tcol(maxlien) + integer(kind=selected_int_kind(16)),INTENT(INOUT),dimension(tnrows)::tcountcol + integer(kind=selected_int_kind(16)),INTENT(INOUT)::tl1(maxlien),tl2(maxlien),tktyp(maxlien) + real*8,INTENT(INOUT)::txtt(maxlien),txjjz(maxlien),txjjxy(maxlien) + + nrows = tnrows + provide nrows + do i=1,maxlien + l1(i)=tl1(i) + l2(i)=tl2(i) + ktyp(i)=tktyp(i) + xtt(i) = txtt(i) + xjjxy(i) = txjjxy(i) + xjjz (i) = txjjz (i) + enddo + ntrou = tntrou + isz = tisz + tcol=0 + tval=0d0 + provide l1 l2 ktyp xtt xjjxy xjjz ntrou +!print *,"l1" +!print *,l1 +!print *,"xjjz" +!print *,xjjz + call unit(tistart, tcountcol,tcol,tval) + + end + + diff --git a/ylogic.irp.f b/ylogic.irp.f new file mode 100644 index 0000000..01356a1 --- /dev/null +++ b/ylogic.irp.f @@ -0,0 +1,94 @@ +! SUBROUTINE ylogic(ideter,yalt,ytrou,yrep1) +BEGIN_PROVIDER [logical, yalt,(maxlien)] +&BEGIN_PROVIDER [logical, ytrou,(maxlien)] +&BEGIN_PROVIDER [logical, yrep1,(maxlien)] + + implicit none + BEGIN_DOC + ! provides the important logical units + ! for the presence of holes (ytrou) and + ! an alternance (yalt) or a rep (yrep) + END_DOC + + integer :: i,il + integer :: ik1,ik2 + integer :: ntr,nalt,naltmax + integer :: itypl(maxlien) + logical :: yw + +! write(6,*)'in elem_diag' + +!---mise a zero + nalt=0 + ntr=0 + yw=.FALSE. + do il=1,nlientot + itypl(il)=0 + yalt(il)=.false. + ytrou(il)=.false. + yrep1(il)=.false. + enddo +!--- + + if(yw)write(6,*)'ideter',(deter(il),il=1,natom) + + do il=1,nlientot + + ik1=iliatom1(il) + ik2=iliatom2(il) + + if(yw)write(6,*)'ik1',ik1,'ik2',ik2 + if(deter(ik1).eq.deter(ik2).and. & + deter(ik1).eq.3)then + yrep1(il)=.true. + endif + + if(deter(ik1).ne.deter(ik2))then + if(deter(ik1).eq.1.and.deter(ik2).eq.2)then + yalt(il)=.true. + nalt=nalt+1 + itypl(il)=1 +! goto 32 + CYCLE + endif + if(deter(ik1).eq.2.and.deter(ik2).eq.1)then + yalt(il)=.true. + nalt=nalt+1 + itypl(il)=2 + goto 32 + CYCLE + endif + 32 continue + if(deter(ik1).eq.1.and.deter(ik2).eq.3)then + ytrou(il)=.true. + ntr=ntr+1 + itypl(il)=3 +! goto 33 + CYCLE + endif + if(deter(ik1).eq.3.and.deter(ik2).eq.1)then + ytrou(il)=.true. + ntr=ntr+1 + itypl(il)=4 +! goto 33 + CYCLE + endif + if(deter(ik1).eq.2.and.deter(ik2).eq.3)then + ytrou(il)=.true. + ntr=ntr+1 + itypl(il)=5 +! goto 33 + CYCLE + endif + if(deter(ik1).eq.3.and.deter(ik2).eq.2)then + ytrou(il)=.true. + ntr=ntr+1 + itypl(il)=6 +! goto 33 + CYCLE + endif + endif +!33 continue + enddo + +END_PROVIDER