Community
cancel
Showing results for
Did you mean:
Highlighted
New Contributor II
11 Views

## Array Sizes

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')

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
```

Tags (1)
5 Replies
Highlighted
Valued Contributor III
11 Views

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 `

Highlighted
Black Belt
11 Views
```    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

Highlighted
New Contributor II
11 Views
```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```

Highlighted
New Contributor II
11 Views

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 ?

Highlighted
New Contributor II
11 Views

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.