-- FreeFem++ v4.14 (mer. 06 mars 2024 16:59:04 CET - git v4.14-1-g2b2052ae) file : svd_test.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // svd_test.edp 2 : // 3 : // Discussion: 4 : // 5 : // Demonstrate the interface between FreeFem++ and LAPACK by calling DGESVD. 6 : // 7 : // Licensing: 8 : // 9 : // This code is distributed under the MIT license. 10 : // 11 : // Modified: 12 : // 13 : // 21 July 2015 14 : // 15 : // Author: 16 : // 17 : // John Burkardt 18 : // 19 : cout << "\n"; 20 : cout << "svd_test():\n"; 21 : cout << " FreeFem++ version\n"; 22 : cout << " Demonstrate the interface to LAPACK by accessin ... : g the DGESVD\n"; 23 : cout << " routine which computes the singular value decom ... : position (SVD).\n"; 24 : 25 : load "lapack" Add lapack interface ... 26 : // 27 : // Here's a function which we don't actually use in this example. 28 : // But it demonstrates how one can take a sparse matrix "Asparse" 29 : // as stored in FreeFem++ and copy it into a dense matrix "Full" 30 : // suitable for working with LAPACK. 31 : // 32 : func real[int,int] arrayFromsparse ( matrix & Asparse ) 33 : { 34 : real[int,int] Full(Asparse.n,Asparse.m); 35 : Full = 0.0; 36 : int[int] I(1); 37 : int[int] J(1); 38 : real[int] C(1); 39 : [I,J,C] = Asparse; 40 : for ( int i = 0; i < I.n; ++i ) 41 : { 42 : Full ( I(i), J(i) ) = C(i); 43 : } 44 : return Full; 45 : } 46 : // 47 : // Here's a function we will need. It wants to turn the vector sigma into an array. 48 : // The dimensions of this array are determined by the column size of U and the row 49 : // size of V. 50 : // 51 : func real[int,int] arrayFromRectDiag ( real[int,int] & U, 52 : real[int] & sigma ,real[int,int] & V ) 53 : { 54 : real[int,int] res(U.n,V.m); 55 : res = 0.0; 56 : for ( int i = 0; i < sigma.n; ++i ) 57 : { 58 : res(i,i) = sigma(i); 59 : } 60 : return res; 61 : } 62 : // 63 : // Here is a small dense example matrix C for which we might want the 64 : // singular value decomposition. 65 : // 66 : real[int,int] C = [ 67 : [ 5.0e-16, -1.863389981e-16, 0.894427191, 0.4472135955 ], 68 : [ 2.916666667e-16, 1.490711985e-16, -5.123475383e-16, 0.1224361692 ], 69 : [ -0.0, 0.0, 0.0, 0.0 ] ]; 70 : 71 : cout << "\n"; 72 : cout << " The C array:" << C << endl; 73 : // 74 : // Declare arrays for SVD decomposition C = U * SIGMA * V' 75 : // 76 : real[int,int] U(1,1); 77 : real[int,int] V(1,1); 78 : real[int] sigma; 79 : // 80 : // Call the LAPACK routine DGESVD for the SVD decomposition. 81 : // "dgesdd" is a special Freefem++-to-LAPACK interface. 82 : // 83 : dgesdd ( C, U, sigma, V ); 84 : // 85 : // The most interesting item is the vector SIGMA of singular values. 86 : // 87 : cout << " Singular value vector sigma:" << sigma << endl; 88 : // 89 : // By the way, the call to dgesdd overwrote the original information in C. 90 : // 91 : cout << "\n"; 92 : cout << "Let's check that original C == U * sigma * V'" << endl; 93 : // 94 : // Call "arrayFromRectDiag" to build a matrix from the vector sigma. 95 : // The arrays U and V are only needed for dimension information. 96 : // 97 : real[int,int] S = arrayFromRectDiag ( U, sigma, V ); 98 : // 99 : // Postmultiply U by SIGMA. 100 : // 101 : real[int,int] US = U * S; 102 : // 103 : // Transpose V. 104 : // 105 : real[int,int] VT = V'; 106 : // 107 : // Postmultiply U*SIGMA by V' 108 : // 109 : real[int,int] USVT = US * VT; 110 : // 111 : // We should have recovered C. 112 : // 113 : cout << "\n"; 114 : cout << " U " << U << endl; 115 : cout << " V' " << VT << endl; 116 : cout << " S = " << S << endl; 117 : cout << " U * S * V' = " << USVT << endl; 118 : // 119 : // Terminate. 120 : // 121 : cout << "\n"; 122 : cout << "svd_test():\n"; 123 : cout << " Normal end of execution.\n"; 124 : 125 : exit ( 0 ); 126 : sizestack + 1024 =1952 ( 928 ) svd_test(): FreeFem++ version Demonstrate the interface to LAPACK by accessing the DGESVD routine which computes the singular value decomposition (SVD). The C array:3 4 5e-16 -1.863389981e-16 0.894427191 0.4472135955 2.916666667e-16 1.490711985e-16 -5.123475383e-16 0.1224361692 0 0 0 0 Singular value vector sigma:3 1.001516052 0.109344467 0 Let's check that original C == U * sigma * V' U 3 3 0.9984679669 -0.05533280354 0 0.05533280354 0.9984679669 0 0 0 1 V' 4 4 4.440892099e-16 -1.785293318e-16 0.8917050274 0.4526169949 2.339274607e-15 1.333395043e-15 -0.4526169949 0.8917050274 -0.894427191 0.4472135955 -2.220446049e-16 1.609823386e-15 -0.4472135955 -0.894427191 -8.881784197e-16 2.053912596e-15 S = 3 4 1.001516052 0 0 0 0 0.109344467 0 0 0 0 0 0 U * S * V' = 3 4 4.299276843e-16 -1.86593552e-16 0.894427191 0.4472135955 2.800048158e-16 1.35682496e-16 -5.471708806e-16 0.1224361692 0 0 0 0 svd_test(): Normal end of execution. current line = 125 exit(0) err code 0 , mpirank 0 CodeAlloc : nb ptr 4058, size :529120 mpirank: 0 Ok: Normal End