Re: Behaviour of svd_lapack

Not sure what to say. To my knowledge, NCL is passing arguments
*directly* to the LAPACK routine.

The 1st two examples at

http://www.ncl.ucar.edu/Document/Functions/Built-in/svd_lapack.shtml

match the IMSL and NAG results.

 From the above documentation:
"Note that the source LAPACK routine (dgesvd) returns V**T, not V.
However, setting optv=1 will result in the return of V."
Maybe this is an issue?

----
I am very busy right now and won't look at this till next week.
Lunde, Bruce N CIV NAVOCEANO, NP1 wrote:
> Hello,
> 
> I have some questions about the behaviour
> or implementation of the NCL routine
> 
> svd_lapack
> 
> I have been testing it against the SVD
> results of MATLAB and Octave (a free MATLAB
> clone).  The results of MATLAB and Octave
> agree (not surprising), but NCL is not
> "quite" the same.
> 
> Bottom line is that svd_lapack
> (0) needs the inputs (U,V) to be float or double,
> (1) seems to not like nrow.GT.ncol ,
> (2) needs the input matrix transposed to get
>       "close to" the same results that MATLAB does,
> (3) does not return the full U matrix,
> (4) returns a transposed U matrix but correct V
>       matrix when compared to MATLAB (after
>       tranposing the input matrix), but
> (5) the numbers it does return DO MATCH those of
>       MATLAB's svd.
> 
> Some of this behaviour may be expected, some not.
> Results of NCL vs MATLAB are at bottom of this
> message.
> 
> Thanks, Bruce
> 
> .............................................
> THE GORY DETAIL ...
> 
> The MATLAB "svd" routine uses the following
> LAPACK routines:
> 
> DGESVD or ZGESVD or SGESVD or CGESVD
> 
> The test case I will go over can be found
> on the MATLAB page at
> 
> http://www.mathworks.com/access/helpdesk/help/techdoc/index.html
> 
> and clicking on the left-side menu item
> "Function Reference".
> 
> **********************
> *** MATLAB RESULTS ***
> **********************
> 
> For the test case, MATLAB creates a matrix such as
> 
> X=[1 2;3 4;5 6;7 8]
> 
> X =
>      1    2
>      3    4
>      5    6
>      7    8
> 
> % size(X) gives 4 2
> 
> As they say on the web page, the statement
> 
> [U,S,V] = svd(X)
> 
> produces
> 
> U =
>     -0.1525   -0.8226   -0.3945   -0.3800
>     -0.3499   -0.4214    0.2428    0.8007
>     -0.5474   -0.0201    0.6979   -0.4614
>     -0.7448    0.3812   -0.5462    0.0407
> 
> S =
>      14.2691         0
>            0    0.6268
>            0         0
>            0         0
> 
> V =
>     -0.6414     0.7672
>     -0.7672    -0.6414
> 
> This can be confirmed in either MATLAB
> or Octave.
> 
> *******************
> *** NCL RESULTS ***
> *******************
> 
> Now for the NCL results (bear with me ... you
> can cut/paste the code below to test it):
> 
> ;-- START NCL CODE --
> 
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> 
> a = (/ (/1,2/), (/3,4/), (/5,6/), (/7,8/) /)*1.0d0
> 
> ; print(dimsizes(a)) gives 4 2 ... matching MATLAB
> 
> dima = dimsizes(a)
> nra  = dima(0)
> nca  = dima(1)
> 
> nrcu = min( (/nra,nca/) )  ; 4, option S
> 
> u = new ( (/nrcu,nca /) , typeof(a) )  ; pre-allocate space 
> v = new ( (/nrcu,nrcu/) , typeof(a) )  ; for returned arrays
> 
> s = svd_lapack (a, "S" , "S", 0,  u, v)
> 
> ;-- END NCL CODE --
> 
> *** OUTPUT MESSAGE ***
> fatal:svd_lapack: The dimensions of 'v' must be
>       min(N,M) x N where (N,M) are the dimensions of 'a'
> fatal:Execute: Error occurred at or near line 16
> 
> ;-- START NCL CODE --
> 
> delete(dima)
> delete(a)
> delete(u)
> delete(v)
> delete(s)
> 
> b = (/ (/1,2/), (/3,4/), (/5,6/), (/7,8/) /)*1.0d0
> a = transpose(b)
> ; print(dimsizes(a)) now gives 2 4
> 
> dima = dimsizes(a)
> nra  = dima(0)
> nca  = dima(1)
> 
> nrcu = min( (/nra,nca/) )  ; 4, option S
> 
> u = new ( (/nrcu,nca /) , typeof(a) )  ; pre-allocate space 
> v = new ( (/nrcu,nrcu/) , typeof(a) )  ; for returned arrays
> 
> s = svd_lapack (a, "S" , "S", 0,  u, v)
> 
> print(s)
> 
> ; (0)     14.26909549926148
> ; (1)     0.6268282324175422
> ;
> ; These numbers MATCH the results of the nonzero
> ; elements of the MATLAB "S" matrix.  It matches
> ; the diagonal of S, and has the singular values
> ; in the correct order.
> 
> write_matrix(u,nca+"F9.4",False)
> 
> ; -0.1525  -0.3499  -0.5474  -0.7448
> ; -0.8226  -0.4214  -0.0201   0.3812
> ;
> ; These numbers MATCH the results of the first
> ; TWO columns of MATLAB U matrix, IF you TRANSPOSE
> ; the NCL u matrix (MATLAB columns 3,4 are missing).
> 
> write_matrix(v,nca+"F9.4",False)
> 
> ; -0.6414   0.7672
> ; -0.7672  -0.6414
> ;
> ; These numbers MATCH the results of the MATLAB
> ; V matrix (WITHOUT the need to transpose).
> 
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk