EVSL  1.1.0 EigenValues Slicing Library
cs_qrsol.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* x=A\b where A can be rectangular; b overwritten with solution */
3 CS_INT cs_qrsol (CS_INT order, const cs *A, CS_ENTRY *b)
4 {
5  CS_ENTRY *x ;
6  css *S ;
7  csn *N ;
8  cs *AT = NULL ;
9  CS_INT k, m, n, ok ;
10  if (!CS_CSC (A) || !b) return (0) ; /* check inputs */
11  n = A->n ;
12  m = A->m ;
13  if (m >= n)
14  {
15  S = cs_sqr (order, A, 1) ; /* ordering and symbolic analysis */
16  N = cs_qr (A, S) ; /* numeric QR factorization */
17  x = cs_calloc (S ? S->m2 : 1, sizeof (CS_ENTRY)) ; /* get workspace */
18  ok = (S && N && x) ;
19  if (ok)
20  {
21  cs_ipvec (S->pinv, b, x, m) ; /* x(0:m-1) = b(p(0:m-1) */
22  for (k = 0 ; k < n ; k++) /* apply Householder refl. to x */
23  {
24  cs_happly (N->L, k, N->B [k], x) ;
25  }
26  cs_usolve (N->U, x) ; /* x = R\x */
27  cs_ipvec (S->q, x, b, n) ; /* b(q(0:n-1)) = x(0:n-1) */
28  }
29  }
30  else
31  {
32  AT = cs_transpose (A, 1) ; /* Ax=b is underdetermined */
33  S = cs_sqr (order, AT, 1) ; /* ordering and symbolic analysis */
34  N = cs_qr (AT, S) ; /* numeric QR factorization of A' */
35  x = cs_calloc (S ? S->m2 : 1, sizeof (CS_ENTRY)) ; /* get workspace */
36  ok = (AT && S && N && x) ;
37  if (ok)
38  {
39  cs_pvec (S->q, b, x, m) ; /* x(q(0:m-1)) = b(0:m-1) */
40  cs_utsolve (N->U, x) ; /* x = R'\x */
41  for (k = m-1 ; k >= 0 ; k--) /* apply Householder refl. to x */
42  {
43  cs_happly (N->L, k, N->B [k], x) ;
44  }
45  cs_pvec (S->pinv, x, b, n) ; /* b(0:n-1) = x(p(0:n-1)) */
46  }
47  }
48  cs_free (x) ;
49  cs_sfree (S) ;
50  cs_nfree (N) ;
51  cs_spfree (AT) ;
52  return (ok) ;
53 }
void * cs_free(void *p)
Definition: cs_malloc.c:22
css * cs_sfree(css *S)
Definition: cs_util.c:54
#define cs
Definition: cs.h:637
#define css
Definition: cs.h:688
CS_INT cs_happly(const cs *V, CS_INT i, double beta, CS_ENTRY *x)
Definition: cs_happly.c:3
#define CS_ENTRY
Definition: cs.h:635
cs * cs_spfree(cs *A)
Definition: cs_util.c:33
CS_INT cs_ipvec(const CS_INT *p, const CS_ENTRY *b, CS_ENTRY *x, CS_INT n)
Definition: cs_ipvec.c:3
CS_INT cs_utsolve(const cs *U, CS_ENTRY *x)
Definition: cs_utsolve.c:3
#define csn
Definition: cs.h:689
#define CS_CSC(A)
Definition: cs.h:659
CS_INT cs_qrsol(CS_INT order, const cs *A, CS_ENTRY *b)
Definition: cs_qrsol.c:3
css * cs_sqr(CS_INT order, const cs *A, CS_INT qr)
Definition: cs_sqr.c:60
void * cs_calloc(CS_INT n, size_t size)
Definition: cs_malloc.c:16
cs * cs_transpose(const cs *A, CS_INT values)
Definition: cs_transpose.c:3
#define CS_INT
Definition: cs.h:627
csn * cs_nfree(csn *N)
Definition: cs_util.c:43
csn * cs_qr(const cs *A, const css *S)
Definition: cs_qr.c:3
CS_INT cs_usolve(const cs *U, CS_ENTRY *x)
Definition: cs_usolve.c:3
CS_INT cs_pvec(const CS_INT *p, const CS_ENTRY *b, CS_ENTRY *x, CS_INT n)
Definition: cs_pvec.c:3