Go to the documentation of this file.00001 #include "cs.h"
00002
00003 int cs_qrsol (int order, const cs *A, double *b)
00004 {
00005 double *x ;
00006 css *S ;
00007 csn *N ;
00008 cs *AT = NULL ;
00009 int k, m, n, ok ;
00010 if (!CS_CSC (A) || !b) return (0) ;
00011 n = A->n ;
00012 m = A->m ;
00013 if (m >= n)
00014 {
00015 S = cs_sqr (order, A, 1) ;
00016 N = cs_qr (A, S) ;
00017 x = cs_calloc (S ? S->m2 : 1, sizeof (double)) ;
00018 ok = (S && N && x) ;
00019 if (ok)
00020 {
00021 cs_ipvec (S->pinv, b, x, m) ;
00022 for (k = 0 ; k < n ; k++)
00023 {
00024 cs_happly (N->L, k, N->B [k], x) ;
00025 }
00026 cs_usolve (N->U, x) ;
00027 cs_ipvec (S->q, x, b, n) ;
00028 }
00029 }
00030 else
00031 {
00032 AT = cs_transpose (A, 1) ;
00033 S = cs_sqr (order, AT, 1) ;
00034 N = cs_qr (AT, S) ;
00035 x = cs_calloc (S ? S->m2 : 1, sizeof (double)) ;
00036 ok = (AT && S && N && x) ;
00037 if (ok)
00038 {
00039 cs_pvec (S->q, b, x, m) ;
00040 cs_utsolve (N->U, x) ;
00041 for (k = m-1 ; k >= 0 ; k--)
00042 {
00043 cs_happly (N->L, k, N->B [k], x) ;
00044 }
00045 cs_pvec (S->pinv, x, b, n) ;
00046 }
00047 }
00048 cs_free (x) ;
00049 cs_sfree (S) ;
00050 cs_nfree (N) ;
00051 cs_spfree (AT) ;
00052 return (ok) ;
00053 }