19 June      2012  11:03:27.504 AM      
 
SVD_DEMO
  FORTRAN90 version
  Demonstrate the singular value decomposition (SVD)
 
  A real MxN matrix A can be factored as:
 
    A = U * S * V'
 
  where
 
    U = MxM orthogonal,
    S = MxN zero except for diagonal,
    V = NxN orthogonal.
 
  The diagonal of S contains only nonnegative numbers
  and these are arranged in descending order.
 
  Matrix row order    M =        5
  Matrix column order N =        8
  Random number SEED    =    123456789
  (Chosen by user.)
 
  We choose a "random" matrix A, with
  integral values between 0 and 10.
 
  The matrix A:
 
  Col         1             2             3             4             5       
  Row
 
    1:   2.00000       1.00000       1.00000       0.00000       9.00000    
    2:   10.0000       3.00000       4.00000       9.00000       8.00000    
    3:   8.00000       1.00000       4.00000       4.00000       1.00000    
    4:   6.00000       0.00000       8.00000       1.00000       0.00000    
    5:   4.00000       6.00000       8.00000       0.00000       3.00000    
 
  Col         6             7             8       
  Row
 
    1:   9.00000       7.00000       6.00000    
    2:   1.00000       6.00000       2.00000    
    3:   4.00000       9.00000       8.00000    
    4:   8.00000       5.00000       4.00000    
    5:   3.00000       9.00000       2.00000    
 
  The orthogonal factor U:
 
  Col         1             2             3             4             5       
  Row
 
    1: -0.421671      0.487547      0.762151     -0.590472E-01 -0.114523E-01
    2: -0.487891     -0.791659      0.229657     -0.328621E-01 -0.285335    
    3: -0.490789     -0.175013E-01 -0.209319      0.530914      0.658131    
    4: -0.406607      0.350797     -0.437254      0.284878     -0.662771    
    5: -0.421845      0.110495     -0.362461     -0.795241      0.214595    
 
  The diagonal factor S:
 
  Col         1             2             3             4             5       
  Row
 
    1:   30.7103       0.00000       0.00000       0.00000       0.00000    
    2:   0.00000       11.2457       0.00000       0.00000       0.00000    
    3:   0.00000       0.00000       9.73453       0.00000       0.00000    
    4:   0.00000       0.00000       0.00000       7.53857       0.00000    
    5:   0.00000       0.00000       0.00000       0.00000       5.17902    
 
  Col         6             7             8       
  Row
 
    1:   0.00000       0.00000       0.00000    
    2:   0.00000       0.00000       0.00000    
    3:   0.00000       0.00000       0.00000    
    4:   0.00000       0.00000       0.00000    
    5:   0.00000       0.00000       0.00000    
 
  The orthogonal factor V:
 
  Col         1             2             3             4             5       
  Row
 
    1: -0.448566     -0.403244     -0.197960      0.308931     -0.140848    
    2: -0.159790     -0.110439     -0.958408E-01 -0.583422      0.208195    
    3: -0.357014      0.836975E-01 -0.570569     -0.285166     -0.406579    
    4: -0.220147     -0.608602      0.813995E-01  0.280262     -0.115517    
    5: -0.307861     -0.145066      0.760172     -0.351410     -0.209274    
    6: -0.350517      0.592593      0.171177      0.192698     -0.466163    
    7: -0.525094      0.111491     -0.636183E-01 -0.207607      0.530699    
    8: -0.322440      0.251308      0.907830E-01  0.447875      0.464135    
 
  Col         6             7             8       
  Row
 
    1: -0.508569      0.374980      0.287168    
    2:  0.343320      0.165197      0.651307    
    3: -0.474527E-01 -0.524943     -0.124675    
    4:  0.661737     -0.116457     -0.176534    
    5: -0.269805     -0.255718     -0.290740E-01
    6:  0.331972      0.350015      0.975347E-01
    7:  0.170335E-01  0.275777     -0.553818    
    8:  0.708900E-02 -0.528221      0.359935    
 
  The product U * S * V':
 
  Col         1             2             3             4             5       
  Row
 
    1:   2.00000       1.00000       1.00000      0.412257E-14   9.00000    
    2:   10.0000       3.00000       4.00000       9.00000       8.00000    
    3:   8.00000       1.00000       4.00000       4.00000       1.00000    
    4:   6.00000      0.388578E-14   8.00000       1.00000      0.621725E-14
    5:   4.00000       6.00000       8.00000      0.249800E-14   3.00000    
 
  Col         6             7             8       
  Row
 
    1:   9.00000       7.00000       6.00000    
    2:   1.00000       6.00000       2.00000    
    3:   4.00000       9.00000       8.00000    
    4:   8.00000       5.00000       4.00000    
    5:   3.00000       9.00000       2.00000    
 
  Frobenius Norm of A, A_NORM =    35.3270    
 
  ABSOLUTE ERROR for A = U*S*V':
  Frobenius norm of difference A-U*S*V' =   0.272896E-13
 
  RELATIVE ERROR for A = U*S*V':
  Ratio of DIF_NORM / A_NORM =   0.772485E-15
 
RANK_ONE_TEST:
  Compare A to the sum of R rank one matrices.
 
         R    Absolute      Relative
              Error         Error
 
         0     35.3270         1.00000    
         1     17.4608        0.494261    
         2     13.3571        0.378100    
         3     9.14616        0.258900    
         4     5.17902        0.146602    
         5    0.272896E-13    0.772485E-15
 
RANK_ONE_PRINT_TEST:
  Print the sums of R rank one matrices.
 
  Rank R =        0
 
  Col         1             2             3             4             5       
  Row
 
    1:   0.00000       0.00000       0.00000       0.00000       0.00000    
    2:   0.00000       0.00000       0.00000       0.00000       0.00000    
    3:   0.00000       0.00000       0.00000       0.00000       0.00000    
    4:   0.00000       0.00000       0.00000       0.00000       0.00000    
    5:   0.00000       0.00000       0.00000       0.00000       0.00000    
 
  Col         6             7             8       
  Row
 
    1:   0.00000       0.00000       0.00000    
    2:   0.00000       0.00000       0.00000    
    3:   0.00000       0.00000       0.00000    
    4:   0.00000       0.00000       0.00000    
    5:   0.00000       0.00000       0.00000    
 
  Rank R =        1
 
  Col         1             2             3             4             5       
  Row
 
    1:   5.80876       2.06922       4.62320       2.85082       3.98668    
    2:   6.72098       2.39418       5.34924       3.29852       4.61276    
    3:   6.76090       2.40840       5.38101       3.31811       4.64016    
    4:   5.60125       1.99530       4.45805       2.74898       3.84426    
    5:   5.81116       2.07008       4.62511       2.85200       3.98833    
 
  Col         6             7             8       
  Row
 
    1:   4.53907       6.79978       4.17549    
    2:   5.25189       7.86762       4.83121    
    3:   5.28309       7.91435       4.85991    
    4:   4.37692       6.55686       4.02632    
    5:   4.54094       6.80258       4.17721    
 
  Rank R =        2
 
  Col         1             2             3             4             5       
  Row
 
    1:   3.59786       1.46371       5.08210     -0.486015       3.19131    
    2:   10.3110       3.37739       4.60410       8.71675       5.90425    
    3:   6.84026       2.43014       5.36454       3.43789       4.66871    
    4:   4.01048       1.55963       4.78823      0.348079       3.27198    
    5:   5.31009       1.93285       4.72911       2.09576       3.80807    
 
  Col         6             7             8       
  Row
 
    1:   7.78814       7.41106       5.55336    
    2: -0.238205E-01   6.87504       2.59388    
    3:   5.16645       7.89241       4.81045    
    4:   6.71467       6.99669       5.01772    
    5:   5.27729       6.94112       4.48948    
 
  Rank R =        3
 
  Col         1             2             3             4             5       
  Row
 
    1:   2.12916      0.752649      0.848949      0.117902       8.83116    
    2:   9.86839       3.16313       3.32853       8.89873       7.60369    
    3:   7.24363       2.62542       6.52714       3.27203       3.11977    
    4:   4.85309       1.96757       7.21683      0.160564E-02  0.363447E-01
    5:   6.00857       2.27101       6.74230       1.80855       1.12589    
 
  Col         6             7             8       
  Row
 
    1:   9.05813       6.93906       6.22689    
    2:  0.358863       6.73281       2.79683    
    3:   4.81766       8.02204       4.62546    
    4:   5.98606       7.26748       4.63130    
    5:   4.67331       7.16559       4.16916    
 
  Rank R =        4
 
  Col         1             2             3             4             5       
  Row
 
    1:   1.99165       1.01235      0.975885     -0.685149E-02   8.98759    
    2:   9.79186       3.30766       3.39918       8.82929       7.69074    
    3:   8.48008      0.290375       5.38581       4.39374       1.71330    
    4:   5.51654      0.714629       6.60442      0.603489     -0.718334    
    5:   4.15654       5.76861       8.45187      0.128384       3.23259    
 
  Col         6             7             8       
  Row
 
    1:   8.97235       7.03148       6.02753    
    2:  0.311125       6.78424       2.68588    
    3:   5.58890       7.19113       6.41801    
    4:   6.39989       6.82163       5.59315    
    5:   3.51809       8.41019       1.48416    
 
  Rank R =        5
 
  Col         1             2             3             4             5       
  Row
 
    1:   2.00000       1.00000       1.00000      0.412257E-14   9.00000    
    2:   10.0000       3.00000       4.00000       9.00000       8.00000    
    3:   8.00000       1.00000       4.00000       4.00000       1.00000    
    4:   6.00000      0.388578E-14   8.00000       1.00000      0.621725E-14
    5:   4.00000       6.00000       8.00000      0.249800E-14   3.00000    
 
  Col         6             7             8       
  Row
 
    1:   9.00000       7.00000       6.00000    
    2:   1.00000       6.00000       2.00000    
    3:   4.00000       9.00000       8.00000    
    4:   8.00000       5.00000       4.00000    
    5:   3.00000       9.00000       2.00000    
 
  Original matrix A:
 
  Col         1             2             3             4             5       
  Row
 
    1:   2.00000       1.00000       1.00000       0.00000       9.00000    
    2:   10.0000       3.00000       4.00000       9.00000       8.00000    
    3:   8.00000       1.00000       4.00000       4.00000       1.00000    
    4:   6.00000       0.00000       8.00000       1.00000       0.00000    
    5:   4.00000       6.00000       8.00000       0.00000       3.00000    
 
  Col         6             7             8       
  Row
 
    1:   9.00000       7.00000       6.00000    
    2:   1.00000       6.00000       2.00000    
    3:   4.00000       9.00000       8.00000    
    4:   8.00000       5.00000       4.00000    
    5:   3.00000       9.00000       2.00000    
 
  The pseudoinverse of A:
 
  Col         1             2             3             4             5       
  Row
 
    1: -0.289305E-01  0.372563E-01  0.159113E-01  0.319513E-01 -0.288546E-01
    2: -0.598832E-02 -0.875057E-03 -0.984532E-02 -0.457147E-01  0.748500E-01
    3: -0.330086E-01  0.996223E-02 -0.539056E-01  0.742211E-01  0.402066E-01
    4: -0.189293E-01  0.534040E-01  0.777347E-02  0.564769E-02 -0.403379E-01
    5:  0.606697E-01  0.460988E-01 -0.625423E-01 -0.210926E-01  0.289757E-02
    6:  0.434277E-01 -0.726661E-02 -0.446685E-01  0.823752E-01 -0.353797E-01
    7:  0.751516E-02 -0.293409E-01  0.624044E-01 -0.624723E-01  0.545672E-01
    8:  0.178959E-01 -0.379505E-01  0.933326E-01 -0.344409E-01 -0.244964E-01
 
PSEUDO_PRODUCT_TEST
 
  The following relations MUST hold:
 
   A  * A+ * A  = A
   A+ * A  * A+ = A+
 ( A  * A+ ) is MxM symmetric;
 ( A+ * A  ) is NxN symmetric
 
  Here are the Frobenius norms of the errors
  in these relationships:
 
   A  * A+ * A  = A              0.203432E-13
   A+ * A  * A+ = A+             0.164835E-15
 ( A  * A+ ) is MxM symmetric;   0.192931E-14
 ( A+ * A  ) is NxN symmetric;   0.210915E-14
 
  In some cases, the matrix A * A+
  may be interesting (if M <= N, then
  it MIGHT look like the identity.)
 
 
  A * A+:
 
  Col         1             2             3             4             5       
  Row
 
    1:   1.00000     -0.277556E-15 -0.333067E-15 -0.166533E-15 -0.166533E-15
    2: -0.513478E-15   1.00000      0.138778E-15 -0.444089E-15  0.631439E-15
    3:  0.138778E-15   0.00000       1.00000     -0.222045E-15  0.832667E-16
    4:  0.416334E-16 -0.222045E-15   0.00000       1.00000      0.319189E-15
    5: -0.763278E-16 -0.124900E-15 -0.194289E-15 -0.555112E-15   1.00000    
 
  In some cases, the matrix A+ * A
  may be interesting (if N <= M, then
  it MIGHT look like the identity.)
 
 
  A+ * A:
 
  Col         1             2             3             4             5       
  Row
 
    1:  0.518282     -0.743781E-01  0.208513      0.430903     -0.329765E-01
    2: -0.743781E-01  0.430641      0.184212     -0.929715E-01  0.153809    
    3:  0.208513      0.184212      0.706639     -0.517415E-01 -0.150665    
    4:  0.430903     -0.929715E-01 -0.517415E-01  0.517377      0.143627    
    5: -0.329765E-01  0.153809     -0.150665      0.143627      0.860968    
    6:  0.957332E-02 -0.235319      0.211651     -0.161698      0.181909    
    7:  0.642907E-01  0.309300      0.765277E-01 -0.769232E-01  0.590150E-01
    8:  0.983159E-01 -0.149601     -0.232075     -0.266541E-02 -0.122698    
 
  Col         6             7             8       
  Row
 
    1:  0.957332E-02  0.642907E-01  0.983159E-01
    2: -0.235319      0.309300     -0.149601    
    3:  0.211651      0.765277E-01 -0.232075    
    4: -0.161698     -0.769232E-01 -0.266541E-02
    5:  0.181909      0.590150E-01 -0.122698    
    6:  0.757771     -0.481642E-01  0.147426    
    7: -0.481642E-01  0.616943      0.344889    
    8:  0.147426      0.344889      0.591378    
 
PSEUDO_LINEAR_SOLVE_TEST
 
  Given:
    b = A * x1
  so that b is in the range of A, solve
    A * x = b
  using the pseudoinverse:
    x2 = A+ * b.
 
  Norm of x1 =    12.9228    
  Norm of x2 =    12.0214    
  Norm of residual =   0.182542E-12
 
  Given:
    b = A' * x1
  so that b is in the range of A', solve
    A' * x = b
  using the pseudoinverse:
    x2 = A+' * b.
 
  Norm of x1 =    13.3791    
  Norm of x2 =    13.3791    
  Norm of residual =   0.127106E-12
 
  For M < N, most systems A'*x=b will not be
  exactly and uniquely solvable, except in the
  least squares sense.
 
  Here is an example:
 
  Given b is NOT in the range of A', solve
    A' * x = b
  using the pseudoinverse:
    x2 = A+' * b.
 
  Norm of x2 =   0.827285E-01
  Norm of residual =   0.857524    
 
COMPARE_LINPACK_LAPACK:
  While the singular values should be identical,
  the orthogonal factors may have some differences.
 
  Frobenius differences:
 
  U(LAPACK) - U(LINPACK):    2.82843    
  S(LAPACK) - S(LINPACK):   0.141273E-13
  V(LAPACK) - V(LINPACK):    3.46695    
 
SVD_DEMO:
  Normal end of execution.
 
19 June      2012  11:03:27.508 AM