Behaviour of svd_lapack

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