• API
• FAQ
• Tools
• Archive
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, &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.

Top