← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/umfpack/umf_blas3_update.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/pysparse/umfpack/umf_blas3_update.c
 
1
/* ========================================================================== */
 
2
/* === UMF_blas3_update ===================================================== */
 
3
/* ========================================================================== */
 
4
 
 
5
/* -------------------------------------------------------------------------- */
 
6
/* UMFPACK Version 4.1 (Apr. 30, 2003), Copyright (c) 2003 by Timothy A.      */
 
7
/* Davis.  All Rights Reserved.  See ../README for License.                   */
 
8
/* email: davis@cise.ufl.edu    CISE Department, Univ. of Florida.            */
 
9
/* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
 
10
/* -------------------------------------------------------------------------- */
 
11
 
 
12
#include "umf_internal.h"
 
13
 
 
14
GLOBAL void UMF_blas3_update
 
15
(
 
16
    WorkType *Work
 
17
)
 
18
{
 
19
    /* ---------------------------------------------------------------------- */
 
20
    /* local variables */
 
21
    /* ---------------------------------------------------------------------- */
 
22
 
 
23
    Entry *L, *U, *C, *LU ;
 
24
    Int k, m, n, d, nb, dc ;
 
25
 
 
26
    DEBUG5 (("In UMF_blas3_update "ID" "ID" "ID"\n",
 
27
        Work->fnpiv, Work->fnrows, Work->fncols)) ;
 
28
 
 
29
    k = Work->fnpiv ;
 
30
    if (k == 0)
 
31
    {
 
32
        /* no work to do */
 
33
        return ;
 
34
    }
 
35
 
 
36
    m = Work->fnrows ;
 
37
    n = Work->fncols ;
 
38
 
 
39
    d = Work->fnr_curr ;
 
40
    dc = Work->fnc_curr ;
 
41
    nb = Work->nb ;
 
42
    ASSERT (d >= 0 && (d % 2) == 1) ;
 
43
    C = Work->Fcblock ;     /* ldc is fnr_curr */
 
44
    L = Work->Flblock ;     /* ldl is fnr_curr */
 
45
    U = Work->Fublock ;     /* ldu is fnc_curr, stored by rows */
 
46
    LU = Work->Flublock ;   /* nb-by-nb */
 
47
 
 
48
#ifndef NDEBUG
 
49
    DEBUG5 (("DO RANK-NB UPDATE of frontal:\n")) ;
 
50
    DEBUG5 (("DGEMM : "ID" "ID" "ID"\n", k, m, n)) ;
 
51
    DEBUG7 (("C  block: ")) ; UMF_dump_dense (C,  d, m, n) ;
 
52
    DEBUG7 (("A  block: ")) ; UMF_dump_dense (L,  d, m, k) ;
 
53
    DEBUG7 (("B' block: ")) ; UMF_dump_dense (U, dc, n, k) ;
 
54
    DEBUG7 (("LU block: ")) ; UMF_dump_dense (LU, nb, k, k) ;
 
55
#endif
 
56
 
 
57
    if (k == 1)
 
58
    {
 
59
 
 
60
#ifdef USE_NO_BLAS
 
61
 
 
62
        /* no BLAS available - use plain C code instead */
 
63
        Int i, j ;
 
64
 
 
65
        /* rank-1 outer product to update the C block */
 
66
        for (j = 0 ; j < n ; j++)
 
67
        {
 
68
            Entry u_j = U [j] ;
 
69
            if (IS_NONZERO (u_j))
 
70
            {
 
71
                Entry *c_ij, *l_is ;
 
72
                c_ij = & C [j*d] ;
 
73
                l_is = & L [0] ;
 
74
#pragma ivdep
 
75
                for (i = 0 ; i < m ; i++)
 
76
                {
 
77
                    /* C [i+j*d]-= L [i] * U [j] */
 
78
                    MULT_SUB (*c_ij, *l_is, u_j) ;
 
79
                    c_ij++ ;
 
80
                    l_is++ ;
 
81
                }
 
82
            }
 
83
        }
 
84
 
 
85
#else
 
86
        BLAS_GER (m, n, L, U, C, d) ;
 
87
 
 
88
#endif /* USE_NO_BLAS */
 
89
 
 
90
    }
 
91
    else
 
92
    {
 
93
 
 
94
#ifdef USE_NO_BLAS
 
95
 
 
96
        /* no BLAS available - use plain C code instead */
 
97
        Int i, j, s ;
 
98
 
 
99
        /* triangular solve to update the U block */
 
100
        for (s = 0 ; s < k ; s++)
 
101
        {
 
102
            for (i = s+1 ; i < k ; i++)
 
103
            {
 
104
                Entry l_is = LU [i+s*nb] ;
 
105
                if (IS_NONZERO (l_is))
 
106
                {
 
107
                    Entry *u_ij, *u_sj ;
 
108
                    u_ij = & U [i*dc] ;
 
109
                    u_sj = & U [s*dc] ;
 
110
#pragma ivdep
 
111
                    for (j = 0 ; j < n ; j++)
 
112
                    {
 
113
                        /* U [i*dc+j] -= LU [i+s*nb] * U [s*dc+j] ; */
 
114
                        MULT_SUB (*u_ij, l_is, *u_sj) ;
 
115
                        u_ij++ ;
 
116
                        u_sj++ ;
 
117
                    }
 
118
                }
 
119
            }
 
120
        }
 
121
 
 
122
        /* rank-k outer product to update the C block */
 
123
        /* C = C - L*U' (U is stored by rows, not columns) */
 
124
        for (s = 0 ; s < k ; s++)
 
125
        {
 
126
            for (j = 0 ; j < n ; j++)
 
127
            {
 
128
                Entry u_sj = U [j+s*dc] ;
 
129
                if (IS_NONZERO (u_sj))
 
130
                {
 
131
                    Entry *c_ij, *l_is ;
 
132
                    c_ij = & C [j*d] ;
 
133
                    l_is = & L [s*d] ;
 
134
#pragma ivdep
 
135
                    for (i = 0 ; i < m ; i++)
 
136
                    {
 
137
                        /* C [i+j*d]-= L [i+s*d] * U [s*dc+j] */
 
138
                        MULT_SUB (*c_ij, *l_is, u_sj) ;
 
139
                        c_ij++ ;
 
140
                        l_is++ ;
 
141
                    }
 
142
                }
 
143
            }
 
144
        }
 
145
 
 
146
#else
 
147
 
 
148
        BLAS_TRSM_RIGHT (n, k, LU, nb, U, dc) ;
 
149
        BLAS_GEMM (m, n, k, L, U, dc, C, d) ;
 
150
 
 
151
#endif /* USE_NO_BLAS */
 
152
 
 
153
    }
 
154
 
 
155
#ifndef NDEBUG
 
156
    DEBUG5 (("RANK-NB UPDATE of frontal done:\n")) ;
 
157
    DEBUG5 (("DGEMM : "ID" "ID" "ID"\n", k, m, n)) ;
 
158
    DEBUG7 (("C  block: ")) ; UMF_dump_dense (C,  d, m, n) ;
 
159
    DEBUG7 (("A  block: ")) ; UMF_dump_dense (L,  d, m, k) ;
 
160
    DEBUG7 (("B' block: ")) ; UMF_dump_dense (U, dc, n, k) ;
 
161
    DEBUG7 (("LU block: ")) ; UMF_dump_dense (LU, nb, k, k) ;
 
162
#endif
 
163
 
 
164
    DEBUG2 (("blas3 "ID" "ID" "ID"\n", k, Work->fnrows, Work->fncols)) ;
 
165
}