commit 79bdcfa221a96adaea0bc7daf56507c9f5f5e465 Author: vijay gopal chilkuri Date: Mon Dec 19 23:26:16 2016 +0100 initial commit 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