//
//Singurar value decomposition. 
//
//(Caution) 
//This is implemented by Neumerical Recipe in 3rd Edition.

//Syntax 
svd variables /
(options)
;

//Options
output:xx     //xx is output of decomposition uwv/rev. (default=rev)


//Example 1
//This example is verification
//handin segment and 3D data.
//last 3 record is predicated by command.
hand x1 - 8/
0.477857956 0.319091616 0.302620266 0.730827651 0.49992968  0.699196259 0.060517367 0.265924156
0.00406606 0.647777417 0.376794333 0.249187712  0.002900597 0.363580356 0.621609026 0.688979167
0.144337631 0.134830656 0.368313458 0.678782104 0.956669842 0.712539768 0.259800963 0.972840097
0.478567329 0.999731925 0.126118713 0.174019096 0.07744955  0.166693911 0.596564893 0.949427215
0.193299312 0.622925328 0.886683209 0.650727717 0.735032509 0.681099043 0.854324683 0.954729952
0.234341936 0.980023763 0.985418455 0.744663123 0.645158682 0.182633857 0.856732619 0.569249277
0.769870409 0.535464261 0.595999283 0.992034483 0.614017449 0.980085681 0.378980325 0.282306509
0.687151085 0.049443131 0.958651386 0.397166528 0.045288575 0.719019813 0.065615683 0.970861944
0.16716339  0.3670079   0.681103811 0.834992462 0.652788948 0.004701586 0.164266618 0.505313569
0.293595641 0.866260829 0.523882553 0.396366979 0.843838125 0.235840098 0.782898702 0.891316259
0.667230538 0.10930395  0.704181293 0.855822605 0.681378311 0.11694119  0.504940983 0.586558526
0.003491564 0.211539023 0.372501881 0.855361624 0.130223932 0.108133626 0.049682264 0.898001752
0.573589174 0.30593098  0.933263556 0.039450163 0.486293053 0.9937177   0.134909845 0.442912248
0.069926217 0.014254536 0.893675212 0.10080844  0.492802744 0.151568796 0.416445236 0.75165572
0.051691148 0.494699089 0.106039949 0.434479561 0.403959694 0.09126256  0.478329934 0.119430798
0.136113097 0.629852067 0.68700719  0.204129824 0.67403004  0.349366535 0.806686839 0.474264635
0.632875461 0.419292717 0.452484941 0.061667694 0.968768146 0.326369333 0.505115927 0.050154877
0.950692755 0.136194283 0.637065192 0.874784215 0.743306093 0.40860952  0.502982031 0.718677059
0.030777893 0.463988385 0.178683319 0.759725941 0.809173151 0.038466912 0.931100804 0.745170621
0.458882408 0.75344032  0.230500757 0.057578732 0.105972171 0.932919597 0.18738989  0.583214801
;

put nfx;

//特異値分解
svd x1-8/
output:uwv
;

get freq@ana;
put vwu;

//対角化
if(type == "w") outrec;
vector w[8]=x1-8;
for(i=1;i<=8;i++) outrec;
for(i=1;i<=8;i++) {
  if(i != #) w[i]=0.0;
}
select w_1-8;
put w;

//v 転置
get vwu;
if(type == "v") outrec;
select x1-8;
transpose;
select col_1-8;
put v;

//w*vT
get w;
mxmult w_1-8 by v;
put wv;

//u*w*vT
get vwu;
if(type == "u") outrec;
mxmult x1-8 by wv;
put svdx;

//検算 => all 0
mxsub x1-8 by nfx;