Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include "matrix_io.h"
- /* C prototype for LAPACK routine DGELS */
- void dgels_(
- const char * trans, /* 'N' or 'T' */
- const int * m, /* rows in A */
- const int * n, /* cols in A */
- const int * nrhs, /* cols in B */
- double * A, /* array A */
- const int * lda, /* leading dimension of A */
- double * B, /* array B */
- const int * ldb, /* leading dimension of B */
- double * work, /* workspace array */
- int * lwork, /* workspace size */
- int * info /* status code */
- );
- /* call_dgels : wrapper for LAPACK's DGELS routine
- Purpose:
- Solves the least-squares problem
- minimize || A*x-b ||_2^2
- using LAPACK's DGELS routine. Upon exit, the input vector b is overwriten, i.e.,
- the leading n elements of b contain the solution x.
- Return value:
- The function returns the output `info` from DGELS with the
- following exceptions: the return value is
- -12 if the input A is NULL and/or the input B is NULL
- -13 if A is not a tall matrix
- -14 if the dimensions of A and b are incompatible
- -15 in case of memory allocation errors.
- */
- int call_dgels(matrix_t * A, vector_t * b){
- /* Input checks */
- if ((A == NULL) || (b == NULL) || (A->A == NULL) || (b->v == NULL)){
- INPUT_ERR;
- return -12;
- }
- if (A->m < A->n){
- fprintf(stderr,"Matrix A is not a tall matrix, so m<n");
- return -13;
- }
- if (A->m != b->n){
- fprintf(stderr, "The dimensions of A and b are incompatible");
- return -14;
- }
- /* Transpose A, because dgels_() uses column major and C uses row major*/
- int m = A->n;
- int n = A->m;
- /* Initialisering and declaration of values */
- char transpose = 'T'; /* for m<n solves the least square problem wanted */
- const int nrhs = 1; /* b has dim(1,m) */
- int lda = (m == 0) ? 1 : m; /* max(1,m) */
- int ldb = (n == 0) ? 1 : n;; /* max(1,m,n) - m>n */
- int lwork = (n == 0 || m ==0) ? 1 : 2*n; /* def. for lwork max(1,2n) (see report) */
- int info; /* Return value */
- /* Workspace allocation */
- double * work = (double*)malloc(lwork*sizeof(double));
- if (work == 0){
- MEM_ERR;
- return -15;
- }
- /* Calling the function */
- dgels_(&transpose, &m, &n, &nrhs, A->A[0], &lda, b->v, &ldb, work, &lwork, &info);
- /* Free the workspace */
- free(work);
- return info;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement