- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 ?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page