-- FreeFem++ v4.6 (Thu Apr 2 15:47:38 CEST 2020 - git v4.6) 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 : // Location: 8 : // 9 : // http://people.sc.fsu.edu/~jburkardt/freefem_src/svd_test/svd_test.edp 10 : // 11 : // Modified: 12 : // 13 : // 21 July 2015 14 : // 15 : // Author: 16 : // 17 : // John Burkardt 18 : // 19 : load "lapack" Add lapack interface ... 20 : // 21 : // Here's a function which we don't actually use in this example. 22 : // But it demonstrates how one can take a sparse matrix "Asparse" 23 : // as stored in FreeFem++ and copy it into a dense matrix "Full" 24 : // suitable for working with LAPACK. 25 : // 26 : func real[int,int] arrayFromsparse ( matrix & Asparse ) 27 : { 28 : real[int,int] Full(Asparse.n,Asparse.m); 29 : Full = 0.0; 30 : int[int] I(1); 31 : int[int] J(1); 32 : real[int] C(1); 33 : [I,J,C] = Asparse; 34 : for ( int i = 0; i < I.n; ++i ) 35 : { 36 : Full ( I(i), J(i) ) = C(i); 37 : } 38 : return Full; 39 : } 40 : // 41 : // Here's a function we will need. It wants to turn the vector sigma into an array. 42 : // The dimensions of this array are determined by the column size of U and the row 43 : // size of V. 44 : // 45 : func real[int,int] arrayFromRectDiag ( real[int,int] & U, 46 : real[int] & sigma ,real[int,int] & V ) 47 : { 48 : real[int,int] res(U.n,V.m); 49 : res = 0.0; 50 : for ( int i = 0; i < sigma.n; ++i ) 51 : { 52 : res(i,i) = sigma(i); 53 : } 54 : return res; 55 : } 56 : cout << "\n"; 57 : cout << "svd_test:\n"; 58 : cout << " Demonstrate the FreeFem++ interface to LAPACK b ... : y accessing\n"; 59 : cout << " the DGESVD routine which computes the singular ... : value decomposition.\n"; 60 : // 61 : // Here is a small dense example matrix C for which we might want the 62 : // singular value decomposition. 63 : // 64 : real[int,int] C = [ 65 : [ 5.0e-16, -1.863389981e-16, 0.894427191, 0.4472135955 ], 66 : [ 2.916666667e-16, 1.490711985e-16, -5.123475383e-16, 0.1224361692 ], 67 : [ -0.0, 0.0, 0.0, 0.0 ] ]; 68 : 69 : cout << "\n"; 70 : cout << " The C array:" << C << endl; 71 : // 72 : // Declare arrays for SVD decomposition C = U * SIGMA * V' 73 : // 74 : real[int,int] U(1,1); 75 : real[int,int] V(1,1); 76 : real[int] sigma; 77 : // 78 : // Call the LAPACK routine DGESVD for the SVD decomposition. 79 : // "dgesdd" is a special Freefem++-to-LAPACK interface. 80 : // 81 : dgesdd ( C, U, sigma, V ); 82 : // 83 : // The most interesting item is the vector SIGMA of singular values. 84 : // 85 : cout << " Singular value vector sigma:" << sigma << endl; 86 : // 87 : // By the way, the call to dgesdd overwrote the original information in C. 88 : // 89 : cout << "\n"; 90 : cout << "Let's check that original C == U * sigma * V'" << endl; 91 : // 92 : // Call "arrayFromRectDiag" to build a matrix from the vector sigma. 93 : // The arrays U and V are only needed for dimension information. 94 : // 95 : real[int,int] S = arrayFromRectDiag ( U, sigma, V ); 96 : // 97 : // Postmultiply U by SIGMA. 98 : // 99 : real[int,int] US = U * S; 100 : // 101 : // Transpose V. 102 : // 103 : real[int,int] VT = V'; 104 : // 105 : // Postmultiply U*SIGMA by V' 106 : // 107 : real[int,int] USVT = US * VT; 108 : // 109 : // We should have recovered C. 110 : // 111 : cout << "\n"; 112 : cout << " U " << U << endl; 113 : cout << " V' " << VT << endl; 114 : cout << " S = " << S << endl; 115 : cout << " U * S * V' = " << USVT << endl; 116 : // 117 : // Terminate. 118 : // 119 : cout << "\n"; 120 : cout << "svd_test:\n"; 121 : cout << " Normal end of execution.\n"; 122 : 123 : sizestack + 1024 =1952 ( 928 ) svd_test: Demonstrate the FreeFem++ interface to LAPACK by accessing the DGESVD routine which computes the singular value decomposition. 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. times: compile 0.005323s, execution 0.000455s, mpirank:0 CodeAlloc : nb ptr 3680, size :480608 mpirank: 0 Ok: Normal End