From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>

Date: Thu, 02 Jul 2009 09:12:51 -0600

Date: Thu, 02 Jul 2009 09:12:51 -0600

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-talkReceived on Thu Jul 02 2009 - 09:12:51 MDT

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