[ncl-talk] Need a solution for calculate Determinant of a (rather) large matrix
Dennis Shea
shea at ucar.edu
Wed Jul 11 22:20:03 MDT 2018
Your det.stub contains:
C NCLFORTSTART
subroutine TSRGT (eps, n, A, it, C, Kp, Lp)
real*8 eps, po, t0
integer n, it
dimension A(n,n), C(n,n), Kp(n), Lp(n)
C NCLEND
C NCLFORTSTART
*FUNCTION* DMGT (eps, n, A)
integer n
real*8 eps, d0
*dimension *A(n,n)
C NCLEND
---------------------------
[1] Comment: A stub file is required **only* *for the fortran
subroutine(s) invoked from NCL.
d = DET::*DMGT*(eps, N, a)
Hence, the TSRGT stub is not needed. No harm is done... just not needed
-----------------------
[2] The errors:
The actual fortran function was typed as *real*8 *[*double precision*] and
had 'A' typed as* real*8 [double precision]*
* real*8 Function* DMGT(eps, n, A)
integer n
* real*8* eps, A(n,n)
or, equivalently
*double precision Function* DMGT(eps, n, A)
integer n
* double precision *eps, A(n,n)
Your stub file had
FUNCTION DMGT (eps, n, A)
*dimension *A(n,n)
When WRAPIT encounters '*dimension*', it will assign fortran default
typing. For a variable starting with 'A' ['a'], that would be *real*4*. You
must be careful about passing consistent argument (variable) types.
C NCLFORTSTART
FUNCTION DMGT (eps, n, A)
integer n
*double precision* eps, d0
*double precision* A(n,n)
C NCLEND
======================
%> WRAPIT det.stub det.f90
%> ncl check.ncl_new
It still did not return the correct answer. :-(
I suggest that you test the 'det.f90' entirely via fortran. I speculate the
problem is with the fortran.
Good luck
On Wed, Jul 11, 2018 at 5:08 AM, Lam Hoang <lamhpvn at gmail.com> wrote:
> Dear Dennis and all,
> I tied to use WRAPIT to wrap a fortran code to compute determinant.
> However, I got a warning of "*implicit declaration of function "*
>
> *tsrgt_" is invalid in C99" when I wrap tsrgt.f90 using command:*
>
>
> WRAPIT det.stub det.f90
>
> It compiled the code still and generated "det.so" file, however the
> determinant result was wrong (det=0.0 compared to det=30. using
> determinant() function) when I run the code in "check.ncl" file
>
> I think there is something wrong with subroutine TSRGT() because when I
> tried to put subroutine TSRGT() in a separated file and retry to wrap and I
> still got the above warning of implicit declaration.
>
> Attachments are all my used files. Could you please have a check and tell
> me what's wrong.
>
> Regards,
> Lam
>
>
> On Wed, Jul 11, 2018 at 9:43 AM, Lam Hoang <lamhpvn at gmail.com> wrote:
>
>> Thank you Dennis for your kind support.
>> I will used WRAPIT and fortran code to compute determinant instead.
>>
>> Thanks again!
>>
>> Lam
>>
>> On Wed, Jul 11, 2018 at 5:34 AM, Dennis Shea <shea at ucar.edu> wrote:
>>
>>> Yes. The 'determinant' function is taking a *very* long time for a
>>> 41x41 square matrix.
>>>
>>> The underlying fortran code uses a partial-pivoting Gaussian elimination
>>> scheme.This method is quick enough for small matrices. Likely, some
>>> convergence criterion is not being met. It has been a long time since I
>>> have looked at this approach. I speculate that it may behave like N^2
>>> (41^2) time,
>>>
>>> For larger matrices, A better method might be a QR decomposition. Still
>>> , this method has issues also.
>>>
>>> The NCL function
>>> https://test.www.ncl.ucar.edu/Document/Functions/Built-in/so
>>> lve_linsys.shtml
>>> uses LU-decomposition to solve a
>>>
>>> The LAPACK numerical analysts have a solid reason for not developing an
>>> *explicit* determinant subroutine. Too many numerical issues.
>>> Often, computing the determinant is *not* what you should be doing to
>>> solve a given problem.
>>>
>>> I will update the 'determinant' documentation to reflect some of the
>>> above issues.
>>>
>>> ===
>>> I suggest calling a fortran or C code that performs an LU decomposition
>>> and compute the determinant if that is what you need.
>>>
>>> Regards
>>>
>>>
>>>
>>> On Tue, Jul 10, 2018 at 8:20 AM, Lam Hoang <lamhpvn at gmail.com> wrote:
>>>
>>>> Please find the attachment for my "x" matrix
>>>> I run determinant(x) and still waiting for the result after 2mins
>>>>
>>>> On Tue, Jul 10, 2018 at 9:13 PM, Lam Hoang <lamhpvn at gmail.com> wrote:
>>>>
>>>>> Dear Dennis,
>>>>> uname -a gives:
>>>>>
>>>>> Darwin MACBooks-MBP 16.7.0 Darwin Kernel Version 16.7.0: Thu Jan 11
>>>>> 22:59:40 PST 2018; root:xnu-3789.73.8~1/RELEASE_X86_64 x86_64
>>>>>
>>>>> and gcc --version gives:
>>>>>
>>>>>
>>>>> Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr
>>>>> --with-gxx-include-dir=/usr/include/c++/4.2.1
>>>>>
>>>>> Apple LLVM version 8.1.0 (clang-802.0.42)
>>>>>
>>>>> Target: x86_64-apple-darwin16.7.0
>>>>>
>>>>> Thread model: posix
>>>>>
>>>>> InstalledDir: /Applications/Xcode.app/Conten
>>>>> ts/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin
>>>>>
>>>>>
>>>>> anything looks suspicious?
>>>>> Regards,
>>>>>
>>>>>
>>>>> On Tue, Jul 10, 2018 at 8:59 PM, Dennis Shea <shea at ucar.edu> wrote:
>>>>>
>>>>>> Normally, NCL would use an LAPACK code to compute a mathematical
>>>>>> quantity. However, LAPACK contains no such subroutine due to issues with
>>>>>> accuracy and stability:
>>>>>>
>>>>>> http://www.netlib.org/lapack/faq.html#_are_there_routines_in
>>>>>> _lapack_to_compute_determinants
>>>>>>
>>>>>> Really, 'they' (applied mathematicians) do not want people
>>>>>> calculating the determinant.
>>>>>> ===
>>>>>>
>>>>>> NCL uses a computer program from:
>>>>>>
>>>>>> An Introduction to Computational Physics, Tao Pang (Cambridge
>>>>>> University Press,1997).
>>>>>> ===
>>>>>>
>>>>>> The use of integer or float values as input should have no effect.
>>>>>> All non-double input is promoted to double prior to input to the fortran
>>>>>> code used to perform the calculations.
>>>>>>
>>>>>> I have no idea why you are encountering extremely long execution
>>>>>> times. My 6.3.0 [MAC] returns a result more or less immediately. What
>>>>>> version of gcc do you have installed? Please send the output from the
>>>>>> following:
>>>>>>
>>>>>> %> uname -a
>>>>>>
>>>>>> %> gcc --version
>>>>>>
>>>>>> ====
>>>>>> That said: I have tried to calculate the determinant with N=41
>>>>>>
>>>>>> N = 41
>>>>>> x = random_uniform(0,5,(/N,N/))
>>>>>> d = determinant(x) ; VERY fast;
>>>>>> instantaneous return value
>>>>>> print(d)
>>>>>>
>>>>>> Variable: d
>>>>>> Type: float
>>>>>> Total Size: 4 bytes
>>>>>> 1 values
>>>>>> Number of Dimensions: 1
>>>>>> Dimensions and sizes: [1]
>>>>>> Coordinates:
>>>>>> (0)
>>>>>>
>>>>>> *1.767279e+31 <==== ???? I have no idea what the correct
>>>>>> answer is*
>>>>>> Ummm, rather large number!!!
>>>>>>
>>>>>> ====
>>>>>> I will open a JIRA trouble ticket and ask the core NCL developers if
>>>>>> everything looks good.
>>>>>>
>>>>>> ====
>>>>>> Again, echoing LAPACK, calculating a determinant has issues.
>>>>>>
>>>>>> Regards
>>>>>>
>>>>>> On Tue, Jul 10, 2018 at 6:34 AM, Lam Hoang <lamhpvn at gmail.com> wrote:
>>>>>>
>>>>>>> Dear all,
>>>>>>> I want to calculate determinant of a 41x41 real (float) matrix. I
>>>>>>> used determinant() function
>>>>>>> https://www.ncl.ucar.edu/Document/Functions/Built-in/determi
>>>>>>> nant.shtml
>>>>>>>
>>>>>>> However it is very slow.
>>>>>>>
>>>>>>> I try to run the sample in the above page, it's quick. Even if I
>>>>>>> replace integer values by float values (i.e replace 1 by 1.0, 0 by 0.0 and
>>>>>>> so on)
>>>>>>> I try to create a 4x4 real matrix, give all members a value of 1.0
>>>>>>> It seems that the determinant() function runs forever and does not
>>>>>>> return any value
>>>>>>>
>>>>>>> In the following illustration, a is sample matrix after replacing
>>>>>>> integer by float values. It takes NCL few seconds to return the value d = 30
>>>>>>>
>>>>>>> x is new 4x4 matrix created by my new() function, it never gives me
>>>>>>> any result.
>>>>>>>
>>>>>>> *Variable: a*
>>>>>>>
>>>>>>> *Type: float*
>>>>>>>
>>>>>>> *Total Size: 64 bytes*
>>>>>>>
>>>>>>> * 16 values*
>>>>>>>
>>>>>>> *Number of Dimensions: 2*
>>>>>>>
>>>>>>> *Dimensions and sizes: [4] x [4]*
>>>>>>>
>>>>>>> *Coordinates: *
>>>>>>>
>>>>>>> *(0) d = 30*
>>>>>>>
>>>>>>>
>>>>>>> *Variable: x*
>>>>>>>
>>>>>>> *Type: float*
>>>>>>>
>>>>>>> *Total Size: 64 bytes*
>>>>>>>
>>>>>>> * 16 values*
>>>>>>>
>>>>>>> *Number of Dimensions: 2*
>>>>>>>
>>>>>>> *Dimensions and sizes: [4] x [4]*
>>>>>>>
>>>>>>> *Coordinates: *
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> Any suggestion or solution for this?
>>>>>>>
>>>>>>> Regards,
>>>>>>>
>>>>>>>
>>>>>>> p/s: I am using NCL v6.3.0
>>>>>>>
>>>>>>> --
>>>>>>> Lam Hoang, PhD
>>>>>>> Chief of Climate Prediction Division,
>>>>>>> Viet Nam National Center of Hydro-Meteorology Forecasting
>>>>>>> National Hydro-Meteorological Services of Vietnam
>>>>>>> 8 Phao Dai Lang street, Lang Thuong, Dong Da, Ha Noi, Vietnam
>>>>>>> <https://maps.google.com/?q=8+Phao+Dai+Lang+street,+Lang+Thuong,+Dong+Da,+Ha+Noi,+Vietnam&entry=gmail&source=g>
>>>>>>> Mob: +84 9 682 34 682
>>>>>>>
>>>>>>> Associate Researcher,
>>>>>>> School of Earth, Atmosphere and Environment
>>>>>>> Monash University, VIC 3800
>>>>>>> Email: lam.hoang at monash.edu
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> ncl-talk mailing list
>>>>>>> ncl-talk at ucar.edu
>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>>
>>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> Lam Hoang, PhD
>>>>> Chief of Climate Prediction Division,
>>>>> Viet Nam National Center of Hydro-Meteorology Forecasting
>>>>> National Hydro-Meteorological Services of Vietnam
>>>>> 8 Phao Dai Lang street, Lang Thuong, Dong Da, Ha Noi, Vietnam
>>>>> <https://maps.google.com/?q=8+Phao+Dai+Lang+street,+Lang+Thuong,+Dong+Da,+Ha+Noi,+Vietnam&entry=gmail&source=g>
>>>>> Mob: +84 9 682 34 682
>>>>>
>>>>> Associate Researcher,
>>>>> School of Earth, Atmosphere and Environment
>>>>> Monash University, VIC 3800
>>>>> Email: lam.hoang at monash.edu
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> Lam Hoang, PhD
>>>> Chief of Climate Prediction Division,
>>>> Viet Nam National Center of Hydro-Meteorology Forecasting
>>>> National Hydro-Meteorological Services of Vietnam
>>>> 8 Phao Dai Lang street, Lang Thuong, Dong Da, Ha Noi, Vietnam
>>>> <https://maps.google.com/?q=8+Phao+Dai+Lang+street,+Lang+Thuong,+Dong+Da,+Ha+Noi,+Vietnam&entry=gmail&source=g>
>>>> Mob: +84 9 682 34 682
>>>>
>>>> Associate Researcher,
>>>> School of Earth, Atmosphere and Environment
>>>> Monash University, VIC 3800
>>>> Email: lam.hoang at monash.edu
>>>>
>>>
>>>
>>
>>
>> --
>> Lam Hoang, PhD
>> Chief of Climate Prediction Division,
>> Viet Nam National Center of Hydro-Meteorology Forecasting
>> National Hydro-Meteorological Services of Vietnam
>> 8 Phao Dai Lang street, Lang Thuong, Dong Da, Ha Noi, Vietnam
>> <https://maps.google.com/?q=8+Phao+Dai+Lang+street,+Lang+Thuong,+Dong+Da,+Ha+Noi,+Vietnam&entry=gmail&source=g>
>> Mob: +84 9 682 34 682
>>
>> Associate Researcher,
>> School of Earth, Atmosphere and Environment
>> Monash University, VIC 3800
>> Email: lam.hoang at monash.edu
>>
>
>
>
> --
> Lam Hoang, PhD
> Chief of Climate Prediction Division,
> Viet Nam National Center of Hydro-Meteorology Forecasting
> National Hydro-Meteorological Services of Vietnam
> 8 Phao Dai Lang street, Lang Thuong, Dong Da, Ha Noi, Vietnam
> <https://maps.google.com/?q=8+Phao+Dai+Lang+street,+Lang+Thuong,+Dong+Da,+Ha+Noi,+Vietnam&entry=gmail&source=g>
> Mob: +84 9 682 34 682
>
> Associate Researcher,
> School of Earth, Atmosphere and Environment
> Monash University, VIC 3800
> Email: lam.hoang at monash.edu
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180711/68ed8f48/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: det.stub
Type: application/octet-stream
Size: 120 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180711/68ed8f48/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: det.f90
Type: application/octet-stream
Size: 4253 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180711/68ed8f48/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: check.ncl_new
Type: application/octet-stream
Size: 1250 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180711/68ed8f48/attachment-0002.obj>
More information about the ncl-talk
mailing list