Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Stability Functions and Pardiso

Valued Contributor III

Dear Guru:

I have been using PARDISO for a while with a structural analysis program from Harrison 1973 Fortran.  Updated of course.

I have no problems with it in the standard solve form. 

I am now looking at high compression members which require a reformulation of the stiffness matrix at each stage using STABILITY FUNCTIONS. These now work and the problems solve in reasonable time. 

I now started a problem - dnl.inp. which is an analysis of a 1 km arch, although that is not relevant. If I just use self weight of the beams it solves nicely, but I started to add the mass of the trucks on the bridge, in this case using a simple method just to develop the model - the 1.1 in the equation increases the deadload. The actual code is in the ELEMENTS.F90 file in the attached solution folder. 

do i = 1,n
        write(*,104)i,nodeloads(i,2),-(nodemass(i)*gr),(-(nodemass(i)*gr) + nodeloads(i,2))
104     Format(' Node : ',I4,' Applied Load :: ',F15.3, ' Self Weight :: ',F15.3, ' Total Load :: ',f15.3)
        nodeloads(i,2) = ZERO - ((nodemass(i)*gr)*1.1)
    end do

If I set 1.1 to 1.0 it runs nicely - about 1.01 it spits out a PARDISO -2 error code at about the 5th iteration.  I cannot see in the help files how to solve for memory size problem. 

Help Please - to run just type dnl at the input. 


0 Kudos
2 Replies
Honored Contributor III

John, you left out the source code of stable.f??, so I had to use the 32-bit OBJ file that you provided. I ran the program with Menu selection 1, then 'dnl' as you wrote.

Every one of the 20 calls to Linsolve() has NNZ not equal to IA(N+1). As long as NNZ > IA(N+1), however, there may be no problem since Pardiso does not directly use NNZ. However, NNZ is used to specify the sizes of some of the arrays, so it should be sufficiently large, and it is worth checking why NNZ is not equal to IA(N+1).

On the 20th call, however, there is a definite problem: N =120, but IA(101:121) = 0. The array IA should have no entry less than the previous one and, in a structure analysis problem, I don't expect any null rows, so each entry should be strictly greater than the previous one. If you do not explicitly ask Pardiso to check the matrix, it won't, but will likely crash because you gave it an impossible task. The error message, of course, is meaningless, and trying to give Pardiso more memory will lead one nowhere.

Beyond that, sorry, I cannot help. Please check the routines where you form the matrix which you then pass to Pardiso.

Please note that in making my diagnosis I completely ignored your description of the "1.1 multiplier". That was not out of disrespect, but simply because I don't understand what that has to do with the zeros in IA. If you are aware of the connection, however, you can probably fix the problem.

0 Kudos
Valued Contributor III

Dear mecej4:

Thanks for the note.  My apologies for stable.f90.  I had developed it in a separate solution folder. I enclose a zip of the folder. 

We have observed that the measured response of the first eigenfrequency has a wider "error" set than would be predicted by linear elasticity.  This becomes really evident on a large English bridge recently.  The geometric change in shape affects the geometry which changes the stiffness matrix. A group from Cambridge LeHigh and Sydney developed the stability function theory, which postulates that the elements in the stiffness matrix can be multiplied by stability functions to allow for the geometric change and the end moments. The alternative formulation is the popular geometric stiffness matrix which is added k(g) matrix added to the stiffness matrix. I prefer the first method, but whilst it is well documented in a few excellent papers - last one attached - I did not find an implementation of this method that I can run on a NUC on site,

You helped me several years ago change Harrison's 73 Fortran code to use the PARDISO solver, a huge help thank you. Harrison had a 2D code sample running the stability functions, and the published functions (attached paper) were tolerably easy to implement in the 3D Harrison Revised Program.  (Balor). 

The S Values are usually around 0.5 to 1, but they can go negative, as that 1.1 factor increases it drives the end moments up which push the S Values negative.  With your great help I now know where to look.  I got fixated on the -2 and should have looked further.  Thanks heaps. 



0 Kudos