SHARE
TWEET

Untitled

a guest Nov 19th, 2019 77 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top