mpskyl.f90

Go to the documentation of this file.

 
 
MODULE mpskyl
    USE mpdef
    IMPLICIT NONE
 
    INTEGER(mpi) :: npar
    INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: lastrowincol
    INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: pargroupstoregions
    INTEGER(mpi), DIMENSION(:), ALLOCATABLE :: regionoffsets
    INTEGER(mpi) :: nthrd
 
END MODULE mpskyl
 
SUBROUTINE mpslan(n,lrows,npgrp,nsparr)
    USE mpskyl
    USE mpdalc
 
    IMPLICIT NONE
    LOGICAL :: lfirst
    INTEGER(mpi) :: i
    INTEGER(mpi) :: ichunk
    INTEGER(mpi) :: j
    INTEGER(mpi) :: ncol
    INTEGER(mpi) :: nrgn
    INTEGER(mpl) :: length
 
    INTEGER(mpi), INTENT(IN) :: n
    INTEGER(mpi), DIMENSION(:), INTENT(IN) :: lrows
    INTEGER(mpi), DIMENSION(:), INTENT(IN) :: npgrp
    INTEGER(mpl), DIMENSION(:,:), INTENT(OUT) :: nsparr
    !$    INTEGER(mpi) :: OMP_GET_MAX_THREADS
    
    npar=n
    print *, ' MPSLAN ', npar
    nthrd=1
    !$ NTHRD=OMP_GET_MAX_THREADS()
    
    length=n
    CALL mpalloc(lastrowincol,length,'last (non-zero) row in column')
    lastrowincol=lrows(:n)
    CALL mpalloc(pargroupstoregions,length,'mapping of parameter groups to regions (in row)')
    length=n+1
    CALL mpalloc(regionoffsets,length,'offsets for parameter group regions (in row)')
    
    ichunk=min((npar+nthrd-1)/nthrd/32+1,256)
    ! analyse sparsity
    ! parallelize row loop
    !$OMP  PARALLEL DO &
    !$OMP  PRIVATE(I,J,NCOL,NRGN,LFIRST) &
    !$OMP  SCHEDULE(DYNAMIC,ICHUNK)
    DO i=1,npar !row
        ncol=0 ! number of columns
        nrgn=0 ! number of regions
        lfirst=.true. 
        DO j=1,i ! loop until diagonal element
            ! inside skyline?
            IF (i <= lastrowincol(j)) THEN
                 IF (lfirst) THEN
                     nrgn=nrgn+1
                     lfirst=.false.
                 END IF
                 ncol=ncol+npgrp(j+1)-npgrp(j)
            ELSE
                lfirst=.true.   
            END IF
        END DO
        nsparr(1,i+1)=2*(nrgn+1)
        nsparr(2,i+1)=ncol*(npgrp(i+1)-npgrp(i))
    END DO 
    !$OMP END PARALLEL DO
        
    !sum up
    nsparr(1,1)=0
    nsparr(2,1)=npgrp(npar+1)-npgrp(1) ! number of diagonal elements
    DO i=1,npar
        nsparr(:,i+1)=nsparr(:,i+1)+nsparr(:,i)
        !print *, ' row ', i, nsparr(:,i+1)
        !print *, ' skyline ', i, lastRowInCol(i), npgrp(i), npgrp(i+1)
    END DO
 
END SUBROUTINE mpslan
 
SUBROUTINE mpslsp(npgrp,nsparr,nsparc)
    USE mpskyl
 
    IMPLICIT NONE
    LOGICAL :: lfirst
    INTEGER(mpi) :: i
    INTEGER(mpi) :: ichunk
    INTEGER(mpi) :: j
    INTEGER(mpl) :: kl
    INTEGER(mpl) :: ll
    INTEGER(mpl) :: l1
    
    INTEGER(mpi), DIMENSION(:), INTENT(IN) :: npgrp
    INTEGER(mpl), DIMENSION(:,:), INTENT(IN) :: nsparr
    INTEGER(mpi), DIMENSION(:), INTENT(OUT) :: nsparc
    
    ichunk=min((npar+nthrd-1)/nthrd/32+1,256)
    ! analyse sparsity
    ! parallelize row loop
    !$OMP  PARALLEL DO &
    !$OMP  PRIVATE(I,J,KL,LL,L1,LFIRST) &
    !$OMP  SCHEDULE(DYNAMIC,ICHUNK)
    DO i=1,npar !row
        kl=nsparr(1,i)    ! pointer to row in NSPARC
        l1=nsparr(2,i)    ! pointer to row in sparse matrix
        ll=l1
        lfirst=.true. 
        DO j=1,i ! loop until diagonal element
            ! inside skyline?
            IF (i <= lastrowincol(j)) THEN
                 IF (lfirst) THEN
                     kl=kl+1
                     nsparc(kl)=int(ll-l1,mpi)
                     kl=kl+1
                     nsparc(kl)=j
                     lfirst=.false.
                 END IF
                 ll=ll+npgrp(j+1)-npgrp(j)
            ELSE
                lfirst=.true.   
            END IF
        END DO
        ! next offset, end marker
        kl=kl+1
        nsparc(kl)=int(ll-l1,mpi)
        kl=kl+1
        nsparc(kl)=npar+1
    END DO 
    !$OMP END PARALLEL DO
    
    WRITE(*,*) ' '
    WRITE(*,*) 'MPSLSP: skyline structure constructed ',nsparr(1,npar+1), ' words'
    WRITE(*,*) 'MPSLSP: dimension parameter of matrix',nsparr(2,1)
    WRITE(*,*) 'MPSLSP: index of last used location',nsparr(2,npar+1)
 
END SUBROUTINE mpslsp
 
SUBROUTINE mpsldc(npgrp,nsparr,nsparc,gmat,nsize,nrank,mon)
    USE mpskyl
    USE mpdalc
 
    IMPLICIT NONE
    INTEGER(mpi) :: i
    INTEGER(mpi) :: ia
    INTEGER(mpi) :: ib
    INTEGER(mpi) :: ik
    INTEGER(mpi) :: ikj
    INTEGER(mpi) :: in
    INTEGER(mpi) :: in1
    INTEGER(mpi) :: ipg
    INTEGER(mpi) :: ir
    INTEGER(mpi) :: j
    INTEGER(mpi) :: jpg
    INTEGER(mpi) :: ja
    INTEGER(mpi) :: jai
    INTEGER(mpi) :: jaj
    INTEGER(mpi) :: jb
    INTEGER(mpi) :: jbi
    INTEGER(mpi) :: jbj
    INTEGER(mpi) :: jk
    INTEGER(mpi) :: jn
    INTEGER(mpi) :: jn1
    INTEGER(mpi) :: jni
    INTEGER(mpi) :: jnj
    INTEGER(mpi) :: jr
    INTEGER(mpi) :: lj
    INTEGER(mpi) :: nrgn
    INTEGER(mpi) :: mpg
    INTEGER(mpl) :: ii
    INTEGER(mpl) :: iij
    INTEGER(mpl) :: is1
    INTEGER(mpl) :: is2
    INTEGER(mpl) :: ki
    INTEGER(mpl) :: kj
    INTEGER(mpl) :: jj
    INTEGER(mpl) :: js1
    INTEGER(mpl) :: js2
    REAL(mpd) :: ratio
    
    INTEGER(mpi), DIMENSION(:), INTENT(IN) :: npgrp
    INTEGER(mpl), DIMENSION(:,:), INTENT(IN) :: nsparr
    INTEGER(mpi), DIMENSION(:), INTENT(IN) :: nsparc
    REAL(mpd), DIMENSION(:), INTENT(INOUT) :: gmat
    INTEGER(mpi), INTENT(IN) :: nsize
    INTEGER(mpi), INTENT(OUT) :: nrank
    INTEGER(mpi), INTENT(IN) :: mon
    
    !$POMP INST BEGIN(mpsldc)
    ! max. used par. group
    mpg=npar
    DO WHILE (npgrp(mpg) > nsize)
        mpg=mpg-1
    END DO
    !print *, ' MPSLDC ', npar, npgrp(npar+1), npgrp(1), nsize, mpg
    !DO ipg=1,npar
    !    print *, ' group ', ipg, npgrp(ipg), npgrp(ipg+1)-npgrp(ipg)
    !END DO
    
    ! parameter groups
    DO ipg=mpg,1,-1
        ! monitoring ?
        IF(mon>0) CALL monpgs(mpg+1-ipg)
        is1=nsparr(1,ipg)
        is2=nsparr(2,ipg)
        in1=int(nsparr(1,ipg+1)-is1,mpi)
        ! inside parameter group
        ia=npgrp(ipg)     ! first (global) row
        ib=npgrp(ipg+1)-1 ! last (global) row
        in=ib-ia+1        ! number of rows 
        ! get mapping par -> region (row-wise expansion)
        nrgn=0
        pargroupstoregions=0
        DO ik=2,in1-2,2
            jpg=nsparc(is1+ik)       ! first group 
            ja=npgrp(jpg)            ! first (global) column
            lj=nsparc(is1+ik-1)      ! region offset 
            jn=nsparc(is1+ik+1)-lj   ! number of columns
            jb=ja+jn-1               ! last (global) column
            nrgn=nrgn+1
            regionoffsets(nrgn+1)=jn
            ! mapping
            DO WHILE (.true.)
                pargroupstoregions(jpg)=nrgn 
                IF (jb < npgrp(jpg+1)) EXIT
                jpg=jpg+1                
            END DO  
        END DO
        !print *, ' listParGroups ', npg, ':', listParGroups(:,:npg) 
        ! sumup
        regionoffsets(1)=0
        DO ir=1,nrgn
            regionoffsets(ir+1)=regionoffsets(ir+1)+regionoffsets(ir)
        END DO 
        !print *, ' region offsets ', nrgn, ':', regionOffsets(:nrgn)     
        ! rows in group 
        DO i=min(ib,nsize),ia,-1
            ! diagonal element
            IF (gmat(i) > 0.0_mpd) THEN
                ! update rank, min, max eigenvalue
                nrank=nrank+1
                !IF (nrank == 1) THEN
                !    evmax=gmat(i)
                !    evmin=gmat(i)
                !ELSE
                !    evmax=max(evmax,gmat(i))
                !    evmin=min(evmin,gmat(i))
                !END IF
                gmat(i)=1.0/gmat(i)  ! (I,I) div !
            END IF
            !print *, ' diag ', i, nrank, gmat(i)
            
            ! loop over (column) groups
            !$OMP  PARALLEL DO &
            !$OMP   DEFAULT(PRIVATE) &
            !$OMP   SHARED(parGroupsToRegions,regionOffsets,npgrp,nsparr,nsparc,gmat,IPG,I,IS1,IS2,IA,IN,IN1) &
            !$OMP   SCHEDULE(DYNAMIC,2)
            DO jpg=1,ipg
                ! region
                ik=pargroupstoregions(jpg)
                IF (ik == 0) cycle
                ii=is2+regionoffsets(ik)*in
                ik=ik+ik
                jr=npgrp(nsparc(is1+ik)) ! first (global) column
                lj=nsparc(is1+ik-1)      ! region offset 
                jn=nsparc(is1+ik+1)-lj   ! number of columns
                ! group
                js1=nsparr(1,jpg)
                js2=nsparr(2,jpg)
                jn1=int(nsparr(1,jpg+1)-js1,mpi)
                ! loop over columns
                ja=npgrp(jpg)            ! first (global) column
                jb=npgrp(jpg+1)-1        ! last (global) column
                lj=(i-ia)*jn-jr+ja       ! offset for row i
                DO j=ja,jb
                    lj=lj+1
                    IF (j >= i) EXIT
                    ratio=gmat(ii+lj)*gmat(i)
                    IF (ratio == 0.0_mpd) cycle
                    jk=2
                    iij=is2+1
                    jj=js2+1
                    DO ikj=2,ik,2
                        ! region in row i
                        jai=npgrp(nsparc(is1+ikj))
                        IF (jai > j) EXIT
                        jni=nsparc(is1+ikj+1)-nsparc(is1+ikj-1)
                        jbi=min(jai+jni-1,j)
                        DO WHILE (.true.)
                            IF (jk >= jn1) THEN
                                print *, ' ERROR in skyline structure ', i, j, ipg, jpg
                                print *, nsparc(is1+1:is1+in1)
                                print *, nsparc(js1+1:js1+jn1)
                                !RETURN
                            END IF
                            jaj=npgrp(nsparc(js1+jk))
                            jnj=nsparc(js1+jk+1)-nsparc(js1+jk-1)
                            jbj=jaj+jnj-1
                            ! (region i) inside (region j)?
                            IF (jaj <= jai.AND.jbj >= jbi) EXIT
                            jk=jk+2
                            jj=jj+jnj*(npgrp(jpg+1)-npgrp(jpg))
                        END DO
                        !print *, ' inside ', jai, jbi, jaj, jbj
                        ! offsets in groups/regions
                        ki=iij+(i-ia)*jni-jai
                        kj=jj+(j-ja)*jnj-jaj
                        ! off diagonal elements
                        gmat(kj+jai:kj+jbi)=gmat(kj+jai:kj+jbi)-gmat(ki+jai:ki+jbi)*ratio
                        ! diagonal element
                        IF (jbi == j) gmat(j)=gmat(j)-gmat(ki+jbi)*ratio
                        iij=iij+in*jni
                    END DO
                END DO
            END DO 
            !$OMP  END PARALLEL DO
                                  
            ! loop over (column) group regions
            ii=is2+1
            DO ik=2,in1-2,2
                jpg=nsparc(is1+ik)       ! first group 
                ja=npgrp(jpg)            ! first (global) column
                lj=nsparc(is1+ik-1)      ! region offset 
                jn=nsparc(is1+ik+1)-lj   ! number of columns
                jb=min(ja+jn-1,i-1)      ! last (global) column
                ki=ii+(i-ia)*jn-ja
                gmat(ki+ja:ki+jb)=gmat(ki+ja:ki+jb)*gmat(i)
                ii=ii+in*jn 
            END DO            
        END DO
    END DO
    !$POMP INST END(mpsldc) 
       
END SUBROUTINE mpsldc
 
SUBROUTINE mpslsv(gij,gvec,nfit)
    USE mpskyl
 
    IMPLICIT NONE
    INTEGER(mpi) :: i
    INTEGER(mpi) :: k
    REAL(mpd) :: dsum
 
    REAL(mpd) :: gij
    REAL(mpd), DIMENSION(:), INTENT(INOUT) :: gvec
    INTEGER(mpi), INTENT(iN) :: nfit
            
    INTERFACE
        FUNCTION gij(i,j)
            USE mpdef
            INTEGER(mpi), INTENT(in) :: i
            INTEGER(mpi), INTENT(in) :: j
        END FUNCTION gij            
    END INTERFACE
        
    DO i=nfit,1,-1 
        dsum=gvec(i)
        DO k=nfit, i+1, -1
            dsum=dsum-gij(k,i)*gvec(k)          ! (K,I)
        END DO
        gvec(i)=dsum
    END DO 
    DO i=1,nfit
        dsum=gvec(i)*gij(i,i)                   ! (I,I)  
        DO k=1,i-1
            dsum=dsum-gij(i,k)*gvec(k)          ! (I,K)
        END DO
        gvec(i)=dsum
    END DO 
 
END SUBROUTINE mpslsv
    
 
!!
SUBROUTINE mpslsv2(npgrp,nsparr,nsparc,gmat,gvec,nsize)
    USE mpskyl
    USE mpdalc
 
    IMPLICIT NONE
    INTEGER(mpi) :: i
    INTEGER(mpi) :: ia
    INTEGER(mpi) :: ib
    INTEGER(mpi) :: ik
    INTEGER(mpi) :: in
    INTEGER(mpi) :: in1
    INTEGER(mpi) :: ipg
    INTEGER(mpi) :: ir
    INTEGER(mpi) :: jpg
    INTEGER(mpi) :: ja
    INTEGER(mpi) :: jai
    INTEGER(mpi) :: jb
    INTEGER(mpi) :: jbi
    INTEGER(mpi) :: jn
    INTEGER(mpi) :: jr
    INTEGER(mpi) :: js
    INTEGER(mpi) :: lj
    INTEGER(mpi) :: irgn
    INTEGER(mpi) :: mpg
    INTEGER(mpl) :: ii
    INTEGER(mpl) :: is1
    INTEGER(mpl) :: is2
    INTEGER(mpl) :: jj
    INTEGER(mpl) :: js1
    INTEGER(mpl) :: js2
    REAL(mpd) :: dsum
    
    INTEGER(mpi), DIMENSION(:), INTENT(IN) :: npgrp
    INTEGER(mpl), DIMENSION(:,:), INTENT(IN) :: nsparr
    INTEGER(mpi), DIMENSION(:), INTENT(IN) :: nsparc
    REAL(mpd), DIMENSION(:), INTENT(IN) :: gmat
    REAL(mpd), DIMENSION(:), INTENT(INOUT) :: gvec
    INTEGER(mpi), INTENT(IN) :: nsize
 
    !$POMP INST BEGIN(mpslsv2b)     
    ! max. used par. group
    mpg=npar
    DO WHILE (npgrp(mpg) > nsize)
        mpg=mpg-1
    END DO
    print *, ' MPSLSV2 ', nsize, npar, mpg
    
    ! 1. backward substitution
    ! get mapping par -> region (column-wise expansion)
    pargroupstoregions=0
    ! parameter groups
    DO ipg=mpg,1,-1
        ! expand (get last region in row)
        irgn=int(nsparr(1,ipg+1)-nsparr(1,ipg),mpi)/2-1
        pargroupstoregions(ipg)=irgn
        regionoffsets(ipg)=int(nsparr(2,ipg+1)-nsparr(2,ipg),mpi)
        !print *, ' ipg, irgn ', ipg, irgn, regionOffsets(ipg)
        ! inside parameter group
        ia=npgrp(ipg)     ! first (global) row
        ib=npgrp(ipg+1)-1 ! last (global) row
        DO i=min(ib,nsize),ia,-1 
            dsum=gvec(i)
            ! loop over rows
            ! $OMP  PARALLEL DO &
            ! $OMP   DEFAULT(PRIVATE) &
            ! $OMP   SHARED(parGroupsToRegions,regionOffsets,npgrp,nsparr,nsparc,gmat,gvec,NPAR,IPG,I,MPG,NSIZE) &
            ! $OMP   REDUCTION(+:DSUM) &
            ! $OMP   SCHEDULE(DYNAMIC,2)
            DO jpg=mpg,ipg,-1
                ! group
                js1=nsparr(1,jpg)
                js2=nsparr(2,jpg)
                ja=npgrp(jpg)            ! first (global) row
                jb=npgrp(jpg+1)-1        ! last (global) row
                js=max(0,jb-nsize)       ! number of rows to skip at end of group
                in=jb-ja+1               ! number of rows
                ! get region with col 'i'
                irgn=pargroupstoregions(jpg)
                DO WHILE (irgn > 0)
                    ik=2*irgn
                    jr=npgrp(nsparc(js1+ik)) ! first (global) column
                    lj=nsparc(js1+ik-1)      ! region offset 
                    jn=nsparc(js1+ik+1)-lj   ! number of columns
                    IF (jr <= i) EXIT
                    irgn=irgn-1
                    pargroupstoregions(jpg)=irgn
                    regionoffsets(jpg)=regionoffsets(jpg)-jn*in
                END DO
                ! not found
                IF (irgn <= 0) cycle
                jai=max(i+1,ja) ! first row
                jbi=jr+jn-1 ! last col in region
                !print *, ' jpg ', jpg, i, irgn, jr, jn , jai, jbi
                IF (jbi < i) cycle
                ! index of end of region
                jj=js2+regionoffsets(jpg)
                ! index of (jai,i)
                ii=jj-jn*(in+ja-jai)-jr+i+1
                dsum=dsum-dot_product(gmat(ii:jj-js*jn:jn),gvec(jai:jb-js))
            END DO
            ! $OMP  END PARALLEL DO
            gvec(i)=dsum
        END DO
    END DO
    !$POMP INST END(mpslsv2b) 
    
    !$POMP INST BEGIN(mpslsv2f) 
    ! 2. forward substitution
    ! parameter groups
    DO ipg=1,mpg
        is1=nsparr(1,ipg)
        is2=nsparr(2,ipg)
        in1=int(nsparr(1,ipg+1)-is1,mpi)
        ! inside parameter group
        ia=npgrp(ipg)     ! first (global) row
        ib=npgrp(ipg+1)-1 ! last (global) row
        in=ib-ia+1        ! number of rows
        ! get mapping par -> region (row-wise expansion)
        irgn=0
        pargroupstoregions=0
        DO ik=2,in1-2,2
            jpg=nsparc(is1+ik)       ! first group 
            ja=npgrp(jpg)            ! first (global) column
            lj=nsparc(is1+ik-1)      ! region offset 
            jn=nsparc(is1+ik+1)-lj   ! number of columns
            jb=ja+jn-1               ! last (global) column
            irgn=irgn+1
            regionoffsets(irgn+1)=jn
            ! mapping
            DO WHILE (.true.)
                pargroupstoregions(jpg)=irgn 
                IF (jb < npgrp(jpg+1)) EXIT
                jpg=jpg+1                
            END DO  
        END DO
        ! sumup
        regionoffsets(1)=0
        DO ir=1,irgn
            regionoffsets(ir+1)=regionoffsets(ir+1)+regionoffsets(ir)
        END DO 
        DO i=ia,min(ib,nsize) 
            dsum=gvec(i)*gmat(i)
            ! columns of row 'i'
            ! loop over (column) groups
            ! $OMP  PARALLEL DO &
            ! $OMP   DEFAULT(PRIVATE) &
            ! $OMP   SHARED(parGroupsToRegions,regionOffsets,npgrp,nsparc,gmat,gvec,IPG,I,IS1,IS2,IA,IN) &
            ! $OMP   REDUCTION(+:DSUM) &
            ! $OMP   SCHEDULE(DYNAMIC,2)
            DO jpg=1,ipg
                ! region
                irgn=pargroupstoregions(jpg)
                IF (irgn == 0) cycle
                ii=is2+regionoffsets(irgn)*in
                ik=2*irgn
                jr=npgrp(nsparc(is1+ik)) ! first (global) column
                lj=nsparc(is1+ik-1)      ! region offset 
                jn=nsparc(is1+ik+1)-lj   ! number of columns
                ! column range
                ja=npgrp(jpg)            ! first (global) column
                jb=npgrp(jpg+1)-1        ! last (global) column
                jbi=min(jb,i-1)
                ii=ii+(i-ia)*jn-jr+1     ! offset for row 'i'
                dsum=dsum-dot_product(gmat(ii+ja:ii+jbi),gvec(ja:jbi))
            END DO
            ! $OMP  END PARALLEL DO
            gvec(i)=dsum
        END DO
    END DO
    !$POMP INST END(mpslsv2f) 
    
END SUBROUTINE mpslsv2