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
This archive was generated by hypermail 2.2.0 : Thu Jul 02 2009 - 11:39:38 MDT