From: Lunde, Bruce N CIV NAVOCEANO, NP1 <bruce.lunde_at_nyahnyahspammersnyahnyah>

Date: Wed, 1 Jul 2009 18:06:38 -0400

Received on Wed Jul 01 2009 - 16:06:38 MDT

Date: Wed, 1 Jul 2009 18:06:38 -0400

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

- application/x-pkcs7-signature attachment: smime.p7s

*
This archive was generated by hypermail 2.2.0
: Thu Jul 02 2009 - 11:39:38 MDT
*