Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Highlighted
##

JohnNichols

New Contributor II

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-15-2016
09:35 AM

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

5 Replies

Highlighted
##

andrew_4619

Valued Contributor III

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-15-2016
10:24 AM

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

jimdempseyatthecove

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-15-2016
11:18 AM

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

JohnNichols

New Contributor II

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-15-2016
11:40 AM

11 Views

Highlighted
##

John_Campbell

New Contributor II

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-15-2016
03:25 PM

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

JohnNichols

New Contributor II

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

06-17-2016
05:10 PM

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.

For more complete information about compiler optimizations, see our Optimization Notice.