Advertisement
Guest User

Untitled

a guest
Nov 19th, 2019
132
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.38 KB | None | 0 0
  1. #include "matrix_io.h"
  2.  
  3. /* C prototype for LAPACK routine DGELS */
  4. void dgels_(
  5. const char * trans, /* 'N' or 'T' */
  6. const int * m, /* rows in A */
  7. const int * n, /* cols in A */
  8. const int * nrhs, /* cols in B */
  9. double * A, /* array A */
  10. const int * lda, /* leading dimension of A */
  11. double * B, /* array B */
  12. const int * ldb, /* leading dimension of B */
  13. double * work, /* workspace array */
  14. int * lwork, /* workspace size */
  15. int * info /* status code */
  16. );
  17.  
  18. /* call_dgels : wrapper for LAPACK's DGELS routine
  19.  
  20. Purpose:
  21. Solves the least-squares problem
  22.  
  23. minimize || A*x-b ||_2^2
  24.  
  25. using LAPACK's DGELS routine. Upon exit, the input vector b is overwriten, i.e.,
  26. the leading n elements of b contain the solution x.
  27.  
  28. Return value:
  29. The function returns the output `info` from DGELS with the
  30. following exceptions: the return value is
  31.  
  32. -12 if the input A is NULL and/or the input B is NULL
  33. -13 if A is not a tall matrix
  34. -14 if the dimensions of A and b are incompatible
  35. -15 in case of memory allocation errors.
  36. */
  37.  
  38. int call_dgels(matrix_t * A, vector_t * b){
  39.  
  40. /* Input checks */
  41. if ((A == NULL) || (b == NULL) || (A->A == NULL) || (b->v == NULL)){
  42. INPUT_ERR;
  43. return -12;
  44. }
  45. if (A->m < A->n){
  46. fprintf(stderr,"Matrix A is not a tall matrix, so m<n");
  47. return -13;
  48. }
  49. if (A->m != b->n){
  50. fprintf(stderr, "The dimensions of A and b are incompatible");
  51. return -14;
  52. }
  53.  
  54. /* Transpose A, because dgels_() uses column major and C uses row major*/
  55. int m = A->n;
  56. int n = A->m;
  57.  
  58. /* Initialisering and declaration of values */
  59. char transpose = 'T'; /* for m<n solves the least square problem wanted */
  60. const int nrhs = 1; /* b has dim(1,m) */
  61. int lda = (m == 0) ? 1 : m; /* max(1,m) */
  62. int ldb = (n == 0) ? 1 : n;; /* max(1,m,n) - m>n */
  63. int lwork = (n == 0 || m ==0) ? 1 : 2*n; /* def. for lwork max(1,2n) (see report) */
  64. int info; /* Return value */
  65.  
  66. /* Workspace allocation */
  67. double * work = (double*)malloc(lwork*sizeof(double));
  68. if (work == 0){
  69. MEM_ERR;
  70. return -15;
  71. }
  72.  
  73. /* Calling the function */
  74. dgels_(&transpose, &m, &n, &nrhs, A->A[0], &lda, b->v, &ldb, work, &lwork, &info);
  75.  
  76. /* Free the workspace */
  77. free(work);
  78.  
  79. return info;
  80. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement