//
//Singurar value decomposition(SVD) 
//SVD generates 3 matrix
// MX=UGV^t
//
// MX.shpe=(20,8) MX is inputed matrix.
//
// U.shape=(20,8)
// G.shape=(8,8) G is diagona.
// V.shape=(8,8) V^t is torasposed V.
//
//(Caution) 
//This is implemented by Neumerical Recipe in 3rd Edition.

//Syntax 
svd variables by n/
(options)
;

//Options
output:xx     //xx is output of decomposition uwv/rev. (default=rev)
              //uwv is SVD  rev is pseudo reverse matirx


//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;

//singular value decomposition
svd x1-8/
output:uwv
;

//reduction of singular value 8->4
svd x1-8/4
output:uwv
;


//verification of SVD
svd x1-8/
putout: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 Transpose
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;

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

//
//Pseudo reverse matrix for non squre matrix
//
get nfx;
svd x1-8/
output:rev
;
get freq@ana;
put revs;

//verify reverse matrix
get nfx;
mxmult x1-8 by revs;
view;
$col1	col2	col3	col4	col5	col6	col7	col8	col9	col10	col11	col12	col13	col14	col15	col16	col17	col18	col19	col20
0.26	-0.06	0.15	-0.06	0.00	0.01	0.28	-0.02	0.11	-0.02	0.00	0.08	0.07	-0.17	0.07	-0.06	0.03	0.05	-0.02	0.13
-0.06	0.31	-0.08	0.12	0.21	0.11	0.06	0.09	-0.16	-0.01	-0.05	0.04	-0.02	0.03	0.06	0.12	-0.20	-0.06	0.12	0.11
0.15	-0.08	0.63	-0.10	0.12	-0.25	-0.03	-0.04	0.07	0.13	-0.08	0.16	0.07	0.09	-0.02	-0.02	0.03	0.00	0.19	0.10
-0.06	0.12	-0.10	0.63	-0.10	0.02	-0.11	0.08	-0.03	0.26	0.02	0.08	-0.10	-0.10	0.00	-0.04	0.05	0.10	0.03	0.26
0.00	0.21	0.12	-0.10	0.29	0.11	0.09	0.06	-0.08	0.01	0.00	0.02	0.09	0.15	0.06	0.19	-0.07	-0.02	0.18	0.03
0.01	0.11	-0.25	0.02	0.11	0.48	0.11	-0.03	0.27	0.09	0.06	0.07	0.04	0.05	0.14	0.17	0.02	-0.09	-0.03	-0.05
0.28	0.06	-0.03	-0.11	0.09	0.11	0.49	0.03	-0.02	-0.13	0.09	-0.02	0.09	-0.21	0.13	0.01	0.01	0.14	0.02	0.12
-0.02	0.09	-0.04	0.08	0.06	-0.03	0.03	0.48	-0.03	-0.08	0.14	0.13	0.19	0.19	-0.16	-0.05	-0.12	0.19	-0.15	0.08
0.11	-0.16	0.07	-0.03	-0.08	0.27	-0.02	-0.03	0.51	0.14	0.06	0.27	0.02	0.04	0.03	-0.04	0.04	-0.06	-0.11	-0.09
-0.02	-0.01	0.13	0.26	0.01	0.09	-0.13	-0.08	0.14	0.30	-0.01	0.03	0.02	0.07	0.04	0.09	0.20	0.00	0.09	0.09
0.00	-0.05	-0.08	0.02	0.00	0.06	0.09	0.14	0.06	-0.01	0.33	0.03	-0.06	0.07	0.00	-0.01	0.08	0.34	0.11	-0.18
0.08	0.04	0.16	0.08	0.02	0.07	-0.02	0.13	0.27	0.03	0.03	0.41	-0.10	0.02	-0.01	-0.13	-0.28	-0.02	0.02	0.00
0.07	-0.02	0.07	-0.10	0.09	0.04	0.09	0.19	0.02	0.02	-0.06	-0.10	0.39	0.15	-0.08	0.11	0.16	-0.04	-0.22	0.18
-0.17	0.03	0.09	-0.10	0.15	0.05	-0.21	0.19	0.04	0.07	0.07	0.02	0.15	0.38	-0.11	0.14	0.03	0.01	0.01	-0.12
0.07	0.06	-0.02	0.00	0.06	0.14	0.13	-0.16	0.03	0.04	0.00	-0.01	-0.08	-0.11	0.14	0.07	0.03	-0.02	0.13	0.00
-0.06	0.12	-0.02	-0.04	0.19	0.17	0.01	-0.05	-0.04	0.09	-0.01	-0.13	0.11	0.14	0.07	0.23	0.12	-0.06	0.10	0.00
0.03	-0.20	0.03	0.05	-0.07	0.02	0.01	-0.12	0.04	0.20	0.08	-0.28	0.16	0.03	0.03	0.12	0.49	0.13	-0.01	0.03
0.05	-0.06	0.00	0.10	-0.02	-0.09	0.14	0.19	-0.06	0.00	0.34	-0.02	-0.04	0.01	-0.02	-0.06	0.13	0.45	0.13	-0.07
-0.02	0.12	0.19	0.03	0.18	-0.03	0.02	-0.15	-0.11	0.09	0.11	0.02	-0.22	0.01	0.13	0.10	-0.01	0.13	0.43	-0.10
0.13	0.11	0.10	0.26	0.03	-0.05	0.12	0.08	-0.09	0.09	-0.18	0.00	0.18	-0.12	0.00	0.00	0.03	-0.07	-0.10	0.40