Array Sizes

JohnNichols · ‎06-15-2016

I could ask this in MKL - but this is more a straight FORTRAN problem.

Before you call PARDISO - you pack the n by n Stiffness matrix into a set of vectors to take up less space. mecej4 showed me in MAGNI how to pack efficiently - the following code does the packing for taking a stiffness matrix from and turns it into the packed vectors, Pardiso needs the length of the vectors essentially NNZ.

I measure the number of degrees of freedom of the problem which gives you a size, but trying to relate NNZ to the degrees of freedom is interesting -

I was thinking just counting the number of non-zero elements in the stiffness matrix and than set NNZ - but I have to loop thru a large matrix

Is there a simpler method, I was just using 9 times the number of degrees of freedom and that usually works, but it did not on the current problem

    subroutine PardisoSolver(GK,nl,nk,ndf,X)


    implicit none
    include 'mkl_pardiso.fi'

    integer nl,nk,ndf

    REAL (KIND=dp) GK(nl,nk),X(ndf)

    REAL (KIND=dp), ALLOCATABLE :: a(:), b(:), c(:,:), atemp(:)
    integer, allocatable :: ja(:), ia(:), iatemp(:), jatemp(:)
    INTEGER :: nnodes = 5
    integer error, iaN
    integer :: nnz , istat, i, j , count
    integer :: size = 0
    integer :: size1 = 0
    integer :: size2 = 0
    integer :: flag = 0
    REAL (KIND=dp) :: delta = 0.000000000000001
    
    nnz = 90*ndf


    !open(16, file="a.mat", STATUS = 'old')

    ! read(16,*)count
    count = ndf

    !nnodes = count
    iaN = count+1

    ALLOCATE (a(nnz), ja(nnz), ia(iaN), jatemp(nnz), iatemp(iaN), b(count), c(count,count), atemp(nnz), STAT=istat)
    IF (istat.NE.0) THEN
        WRITE (*, *) '*** Could not allocate some arrays in LINSOLVE'
        STOP
    END IF

    b = 1.0D0
    a = 0.0d0
    c = 0.0d0
    iatemp = 0
    jatemp = 0
    atemp = 0
    ia = 0
    ja = 0
    do 200 i = 1, count
        flag = 0
        do 300 j = 1, count
            c(i,j) = gk(i,j)
            if(abs(c(i,j)) .gt. delta) then
                size = size + 1
                if(size .gt. nnz) then
                    Write(*,2030)nnz, ndf
2030                Format( "        NNZ is set at :: ",I6, "    NDF is :: ",I6)                    
                    stop    '        Insufficient memory allocation in PARDISO to solve the inversion problem'
                endif
                jatemp(size) = j
                iatemp(size1 + 1) = size+1
                atemp(size) = c(i,j)
                if(flag .eq. 0) then
                    size1 = size1 + 1
                    iatemp(size1) = size
                    flag = 1
                endif
            endif
300 end do
400 format(5(1x,F10.3))
200 end do

andrew_4619 · ‎06-15-2016

count is not a great name for a variable BTW as it is the name of an instrinsic. use of count could remove a loop but if the slice is not contiguous could involve temp arrays being created....

no_zero_count = count( C(:,j)  >  delta )  !  count values in a col

jimdempseyatthecove · ‎06-15-2016

    integer :: size = 0
    integer :: size1 = 0
    integer :: size2 = 0
    integer :: flag = 0

*** Unlike C/C++, those values are 0 only on first call (not re-initialized).

The above is equivalent to the C/C++

    static int size = 0;
    static int size1 = 0;
    static int size2 = 0;
    static int flag = 0;

See: initialization expressions, type declarations

A variable (or variable subobject) can only be initialized once in an executable program.

What you need to do (if the subroutine is called more than once), is to declare the variables without initialization, then in the code section at start of subroutine, initialize the values.

Jim Dempsey

JohnNichols · ‎06-15-2016

subroutine PardisoSolver(GK,nl,nk,ndf,X)


    implicit none
    include 'mkl_pardiso.fi'

    integer nl,nk,ndf

    REAL (KIND=dp) GK(nl,nk),X(ndf),GKTest(nl,nk)

    REAL (KIND=dp), ALLOCATABLE :: a(:), b(:), c(:,:), atemp(:)
    integer, allocatable :: ja(:), ia(:), iatemp(:), jatemp(:)
    INTEGER :: nnodes = 5
    integer error, iaN
    integer :: nnz , istat, i, j !count
    integer :: size = 0
    integer :: size1 = 0
    integer :: size2 = 0
    integer :: flag = 0
    REAL (KIND=dp) :: delta = 0.000000000000001
    integer no_zero_count(nl), numA
    
    GKTest = 0.0d0
    no_zero_count = count( gk .NE. GKTEST, DIM=1 )  !  count values in a col 
    numA = sum(no_zero_count)
    
    write(*,2040)numA
2040 Format('--------------------------------------------------------------------------------------------------------------',///&
            '           Count of the non-zero elements in the stiffness matrix :: ',i6,//)
       
    nnz =  numA+10     !90*ndf
    iaN = ndf+1

    ALLOCATE (a(nnz), ja(nnz), ia(iaN), jatemp(nnz), iatemp(iaN), b(ndf), c(ndf,ndf), atemp(nnz), STAT=istat)
    IF (istat.NE.0) THEN
        WRITE (*, *) '*** Could not allocate some arrays in LINSOLVE'
        STOP
    END IF

    b = 1.0D0
    a = 0.0d0
    c = 0.0d0
    iatemp = 0
    jatemp = 0

John_Campbell · ‎06-15-2016

John,

If you want to find the number of non zero values in GK, why not count them. I am not familiar with GK but assuming nl and nk must be >= ndf, you could write:

!
    if ( nl < ndf ) write (*,*) 'nl is invalid',nl,ndf    
    if ( nk < ndf ) write (*,*) 'nk is invalid',nk,ndf    
    nnz = 0
    do j = 1, ndf
        do i = 1, ndf
            if ( abs(gk(i,j)) .gt. delta) nnz = nnz + 1
        end do
    end do

Note that the inner loop is do i, for sequential scanning of memory.

You also state "I was thinking just counting the number of non-zero elements in the stiffness matrix and than set NNZ - but I have to loop thru a large matrix". If GK is a large matrix, you appear to be wasting a lot of memory, as you store "GK", "c" which is a copy of GK and the sparse form in atemp, jatemp and iatemp. If ndf is truly large, your approach should be to do a scan of the matrix structure to estimate iatemp and then assemble the matrix directly into atemp and jatemp.

Is GK symmetric ?

JohnNichols · ‎06-17-2016

No GK is weird it also holds the load vectors -- as originally written and I have not fixed it.

The count function worked perfectly - I just need to make sure I have enough spaces in the vectors -- better than guessing

PARDISO as Intel has the sample is pretty limited in terms of sending in different size problems through their sample code, it has just been a matter of working out how to make it so it knew the problem size coming in automatically - now it is perfect.

It takes a while to learn the basics but it is worth it - superfast.