[ncl-talk] Need a solution for calculate Determinant of a (rather) large matrix
Lam Hoang
lamhpvn at gmail.com
Wed Jul 11 05:08:49 MDT 2018
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
> 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
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/e8b21120/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: det.f90
Type: application/octet-stream
Size: 4191 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180711/e8b21120/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: det.stub
Type: application/octet-stream
Size: 250 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180711/e8b21120/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: check.ncl
Type: application/octet-stream
Size: 369 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180711/e8b21120/attachment-0002.obj>
More information about the ncl-talk
mailing list