From: Andrew Mai <mai_at_nyahnyahspammersnyahnyah>

Date: Thu, 02 Jul 2009 14:10:19 -0600

Date: Thu, 02 Jul 2009 14:10:19 -0600

Lunde, Bruce N CIV NAVOCEANO, NP1 wrote:

*> I have some questions about the behaviour
*

*> or implementation of the NCL routine
*

*>
*

*> svd_lapack
*

*>
*

You say you have questions but do not ask any here but rather make a

series of observations. I can only say whether your statements are true

or not.

*> 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,
*

*>
*

True.

*> (1) seems to not like nrow.GT.ncol ,
*

*>
*

True, and I do not see any reason for this. The underlying LAPACK

routines can calculate singular values and singular vectors (left and

right) for any M,N > 0.

*> (2) needs the input matrix transposed to get
*

*> "close to" the same results that MATLAB does,
*

*>
*

This is purely a notational issue. Ordinarily we talk about MxN matrices

which have M rows and N columns. These are referred to in NCL as NxM

arrays. Yeah, I agree, it's quirky if you come from a Fortran / MATLAB

background.

*> (3) does not return the full U matrix,
*

*>
*

Only the first 2 columns of U have meaning, even in the LAPACK and

MATLAB versions of this problem. These are the left singular vectors. So

in this sense, NCL *does* return the full U matrix.

When calling the underlying LAPACK routines, there are two

single-character arguments called JOBU and JOBV. The answer you got with

MATLAB sets these two arguments to "A" while NCL sets them to "S" which

means that only the first min(M,N) columns of U are returned. (NCL does

not implement the "A" option through the user interface.)

*> (4) returns a transposed U matrix but correct V
*

*> matrix when compared to MATLAB (after
*

*> tranposing the input matrix),
*

Again, notation only.

*> (5) but the numbers it does return DO MATCH those of
*

*> MATLAB's svd.
*

*>
*

Well, they'd better match or something is seriously wrong.

*> a = (/ (/1,2/), (/3,4/), (/5,6/), (/7,8/) /)*1.0d0
*

*>
*

*> ; print(dimsizes(a)) gives 4 2 ... matching MATLAB
*

*>
*

But you need 2 4 to match MATLAB.

My version of the NCL program (which, note, does get the right answer

when you interpret rows and columns correctly):

a = new((/2,4/),"double")

a(0,0) = 1.0d0

a(1,0) = 2.0d0

a(0,1) = 3.0d0

a(1,1) = 4.0d0

a(0,2) = 5.0d0

a(1,2) = 6.0d0

a(0,3) = 7.0d0

a(1,3) = 8.0d0

print(dimsizes(a))

dima = dimsizes(a)

nca = dima(0)

nra = dima(1)

nrcu = min( (/nra,nca/) ) ; 4, option S

u = new ( (/nrcu,nra /) , typeof(a) ) ; pre-allocate space

v = new ( (/nrcu,nrcu/) , typeof(a) ) ; for returned arrays

s = svd_lapack (a, "S" , "S", 0, u, v)

More questions, comments or observations are welcome,

Andy

_______________________________________________

ncl-talk mailing list

List instructions, subscriber options, unsubscribe:

http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Thu Jul 02 2009 - 14:10:19 MDT

*
This archive was generated by hypermail 2.2.0
: Tue Jul 07 2009 - 11:13:18 MDT
*