← Back to branch summary

~ubuntu-branches/ubuntu/natty/suitesparse/natty

« back to all changes in this revision

Viewing changes to CXSparse/Source/cs_qr.c

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2006-12-22 10:16:15 UTC
  • Revision ID: james.westby@ubuntu.com-20061222101615-2ohaj8902oix2rnk
Tags: upstream-2.3.1
Import upstream version 2.3.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
CXSparse/Source/cs_qr.c
 
1
#include "cs.h"
 
2
/* sparse QR factorization [V,beta,pinv,R] = qr (A) */
 
3
csn *cs_qr (const cs *A, const css *S)
 
4
{
 
5
    CS_ENTRY *Rx, *Vx, *Ax, *Beta, *x ;
 
6
    CS_INT i, k, p, m, n, vnz, p1, top, m2, len, col, rnz, *s, *leftmost, *Ap, *Ai,
 
7
        *parent, *Rp, *Ri, *Vp, *Vi, *w, *pinv, *q ;
 
8
    cs *R, *V ;
 
9
    csn *N ;
 
10
    if (!CS_CSC (A) || !S) return (NULL) ;
 
11
    m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
 
12
    q = S->q ; parent = S->parent ; pinv = S->pinv ; m2 = S->m2 ;
 
13
    vnz = S->lnz ; rnz = S->unz ; leftmost = S->leftmost ;
 
14
    w = cs_malloc (m2+n, sizeof (CS_INT)) ;         /* get CS_INT workspace */
 
15
    x = cs_malloc (m2, sizeof (CS_ENTRY)) ;         /* get CS_ENTRY workspace */
 
16
    N = cs_calloc (1, sizeof (csn)) ;               /* allocate result */
 
17
    if (!w || !x || !N) return (cs_ndone (N, NULL, w, x, 0)) ;
 
18
    s = w + m2 ;                                    /* s is size n */
 
19
    for (k = 0 ; k < m2 ; k++) x [k] = 0 ;          /* clear workspace x */
 
20
    N->L = V = cs_spalloc (m2, n, vnz, 1, 0) ;      /* allocate result V */
 
21
    N->U = R = cs_spalloc (m2, n, rnz, 1, 0) ;      /* allocate result R */
 
22
    N->B = Beta = cs_malloc (n, sizeof (CS_ENTRY)) ;  /* allocate result Beta */
 
23
    if (!R || !V || !Beta) return (cs_ndone (N, NULL, w, x, 0)) ;
 
24
    Rp = R->p ; Ri = R->i ; Rx = R->x ;
 
25
    Vp = V->p ; Vi = V->i ; Vx = V->x ;
 
26
    for (i = 0 ; i < m2 ; i++) w [i] = -1 ; /* clear w, to mark nodes */
 
27
    rnz = 0 ; vnz = 0 ;
 
28
    for (k = 0 ; k < n ; k++)               /* compute V and R */
 
29
    {
 
30
        Rp [k] = rnz ;                      /* R(:,k) starts here */
 
31
        Vp [k] = p1 = vnz ;                 /* V(:,k) starts here */
 
32
        w [k] = k ;                         /* add V(k,k) to pattern of V */
 
33
        Vi [vnz++] = k ;
 
34
        top = n ;
 
35
        col = q ? q [k] : k ;
 
36
        for (p = Ap [col] ; p < Ap [col+1] ; p++)   /* find R(:,k) pattern */
 
37
        {
 
38
            i = leftmost [Ai [p]] ;         /* i = min(find(A(i,q))) */
 
39
            for (len = 0 ; w [i] != k ; i = parent [i]) /* traverse up to k */
 
40
            {
 
41
                s [len++] = i ;
 
42
                w [i] = k ;
 
43
            }
 
44
            while (len > 0) s [--top] = s [--len] ; /* push path on stack */
 
45
            i = pinv [Ai [p]] ;             /* i = permuted row of A(:,col) */
 
46
            x [i] = Ax [p] ;                /* x (i) = A(:,col) */
 
47
            if (i > k && w [i] < k)         /* pattern of V(:,k) = x (k+1:m) */
 
48
            {
 
49
                Vi [vnz++] = i ;            /* add i to pattern of V(:,k) */
 
50
                w [i] = k ;
 
51
            }
 
52
        }
 
53
        for (p = top ; p < n ; p++) /* for each i in pattern of R(:,k) */
 
54
        {
 
55
            i = s [p] ;                     /* R(i,k) is nonzero */
 
56
            cs_happly (V, i, Beta [i], x) ; /* apply (V(i),Beta(i)) to x */
 
57
            Ri [rnz] = i ;                  /* R(i,k) = x(i) */
 
58
            Rx [rnz++] = x [i] ;
 
59
            x [i] = 0 ;
 
60
            if (parent [i] == k) vnz = cs_scatter (V, i, 0, w, NULL, k, V, vnz);
 
61
        }
 
62
        for (p = p1 ; p < vnz ; p++)        /* gather V(:,k) = x */
 
63
        {
 
64
            Vx [p] = x [Vi [p]] ;
 
65
            x [Vi [p]] = 0 ;
 
66
        }
 
67
        Ri [rnz] = k ;                     /* R(k,k) = norm (x) */
 
68
        Rx [rnz++] = cs_house (Vx+p1, Beta+k, vnz-p1) ; /* [v,beta]=house(x) */
 
69
    }
 
70
    Rp [n] = rnz ;                          /* finalize R */
 
71
    Vp [n] = vnz ;                          /* finalize V */
 
72
    return (cs_ndone (N, NULL, w, x, 1)) ;  /* success */
 
73
}