Advertisement
Guest User

Untitled

a guest
Dec 22nd, 2014
147
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 4.57 KB | None | 0 0
  1. apt-file search clapack.h
  2.  
  3. libatlas-dev: /usr/include/atlas/clapack.h
  4. libfreefem++-dev: /usr/include/freefem++/clapack.h
  5.  
  6. apt-get install libfreefem++-dev
  7.  
  8. apt-cache search lapack
  9.  
  10. liblapack-dev - library of linear algebra routines 3 - static version
  11. liblapack3gf - library of linear algebra routines 3 - shared version
  12.  
  13. #include <freefem++/clapack.h>
  14.  
  15. all: lapack_install lib blas_testing lapack_testing
  16.  
  17. all: lapack_install lib
  18.  
  19. make
  20.  
  21. /**
  22. * svd_demo.cpp
  23. *
  24. * Given that you put version 3.5.0 into /opt/lapack/ compile this with:
  25. * g++ svd_demo.cpp -I"/opt/lapack/lapack-3.5.0/lapacke/include"
  26. * -L"/opt/lapack/lapack-3.5.0" -llapacke -llapack -lblas -lcblas
  27. * The order of included libraries is important!
  28. */
  29.  
  30. #include <iostream>
  31. #include <string>
  32. #include <sstream>
  33. #include <cstdlib>
  34. #include <cblas.h>
  35. #include <lapacke.h>
  36.  
  37. using namespace std;
  38.  
  39. typedef double value;
  40.  
  41. /** Column major style! */
  42. string matrix2string(int m, int n, value* A)
  43. {
  44. ostringstream oss;
  45. for (int j=0;j<m;j++)
  46. {
  47. for (int k=0;k<n;k++)
  48. {
  49. oss << A[j+k*m] << "t";
  50. }
  51. oss << endl;
  52. }
  53. return oss.str();
  54. }
  55.  
  56. int main(int argc, char** argv)
  57. {
  58. //> Part 1. Decomposition. -----------------------------------------
  59. char jobu = 'A'; // Return the complete matrix U
  60. char jobvt = 'A'; // Return the complete matrix VT
  61. int mA = 2;
  62. int nA = 3;
  63. int lda = 2;
  64. int ldu = 2;
  65. int ldvt = 3;
  66. int lwork = 81;
  67. int info = 0;
  68. value* A = (value*)malloc(mA*nA*sizeof(value));
  69. value* U = (value*)malloc(mA*mA*sizeof(value));
  70. value* VT = (value*)malloc(nA*nA*sizeof(value));
  71. value* Svec = (value*)malloc(3*sizeof(value));
  72. value* work = (value*)malloc(lwork*sizeof(value));
  73.  
  74. A[0] = 1; A[2] = 2; A[4] = 4;
  75. A[1] = 0; A[3] = 0; A[5] = 4;
  76.  
  77. cout << "Matrix A (will be overwritten, as is documented):" << endl <<
  78. matrix2string(mA,nA,A);
  79.  
  80. // Citing lapacke.h
  81. //lapack_int LAPACKE_dgesvd(int matrix_order, char jobu, char jobvt,
  82. // lapack_int m, lapack_int n, double* a,
  83. // lapack_int lda, double* s, double* u, lapack_int ldu,
  84. // double* vt, lapack_int ldvt, double* superb);
  85.  
  86. info = LAPACKE_dgesvd(LAPACK_COL_MAJOR, jobu, jobvt, mA, nA, A, lda, Svec, U, ldu, VT, ldvt, work);
  87. cout << "Ran dgesvd. Let's see ..." << endl <<
  88. "U:" << endl << matrix2string(mA,mA,U) <<
  89. "Svec:" << endl << matrix2string(1,nA,Svec) <<
  90. "VT:" << endl << matrix2string(nA,nA,VT) <<
  91. "Info Code: " << info << endl << endl <<
  92. "All is well." << endl;
  93. //< ----------------------------------------------------------------
  94. //> Part 2. Checking the result. -----------------------------------
  95. value* S = (value*)malloc(mA*nA*sizeof(value));
  96. S[0] = Svec[0]; S[2] = 0 ; S[4] = 0 ;
  97. S[1] = 0 ; S[3] = Svec[1]; S[5] = 0 ;
  98.  
  99. // Citing cblas.h
  100. // void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
  101. // const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
  102. // const int K, const double alpha, const double *A,
  103. // const int lda, const double *B, const int ldb,
  104. // const double beta, double *C, const int ldc);
  105.  
  106. // work := S*VT; (2x3)=(2x3)*(3x3)
  107. cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,mA,nA,nA,1,S,lda,VT,ldvt,0,work,lda) ;
  108. cout << "Step 1: S*VT" << endl << matrix2string(2,3,work);
  109.  
  110. // A := U*work; (2x2)*(2x3)
  111. cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,mA,nA,mA,1,U,ldu,work,lda,0,A,lda);
  112. cout << "A := U*S*VT:" << endl << matrix2string(mA,nA,A) << endl;
  113. //< ----------------------------------------------------------------
  114. free(A); free(U); free(VT); free(Svec); free(work); free(S);
  115. return EXIT_SUCCESS;
  116. }
  117.  
  118. 1 2 4
  119. 0 0 4
  120. Ran dgesvd. Let's see ...
  121. U:
  122. -0.759729 -0.65024
  123. -0.65024 0.759729
  124. Svec:
  125. 5.89017 1.51851 0
  126. VT:
  127. -0.128982 -0.257965 -0.957506
  128. -0.42821 -0.856419 0.288414
  129. -0.894427 0.447214 -7.48099e-18
  130. Info Code: 0
  131.  
  132. All is well.
  133. Step 1: S*VT
  134. -0.759729 -1.51946 -5.63988
  135. -0.65024 -1.30048 0.437958
  136. A := U*S*VT:
  137. 1 2 4
  138. -9.63558e-16 -4.86265e-17 4
  139.  
  140. libblas-dev - Basic Linear Algebra Subroutines 3, static library
  141. libblas3gf - Basic Linear Algebra Reference implementations, shared library
  142. libopenblas-dev - Optimized BLAS (linear algebra) library based on GotoBLAS2
  143.  
  144. BLASLIB = /usr/lib/openblas-base/libopenblas.a
  145.  
  146. sudo apt-get install libblas-dev checkinstall
  147. sudo apt-get install libblas-doc checkinstall
  148. sudo apt-get install liblapacke-dev checkinstall
  149. sudo apt-get install liblapack-doc checkinstall
  150.  
  151. g++ svd_demo.cpp -I"/usr/include" -L"/usr/lib" -llapacke -lblas
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement