← Back to branch summary

~ubuntu-branches/ubuntu/gutsy/blender/gutsy-security

« back to all changes in this revision

Viewing changes to intern/opennl/superlu/smyblas2.c

  • Committer: Bazaar Package Importer
  • Author(s): Florian Ernst
  • Date: 2005-11-06 12:40:03 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20051106124003-3pgs7tcg5rox96xg
Tags: 2.37a-1.1
http://bugs.debian.org/333958
* Non-maintainer upload.
* Split out parts of 01_SConstruct_debian.dpatch again: root_build_dir
  really needs to get adjusted before the clean target runs - closes: #333958,
  see #288882 for reference

Show diffs side-by-side

added added

removed removed

Lines of Context:
intern/opennl/superlu/smyblas2.c
 
1
 
 
2
 
 
3
/*
 
4
 * -- SuperLU routine (version 2.0) --
 
5
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 
6
 * and Lawrence Berkeley National Lab.
 
7
 * November 15, 1997
 
8
 *
 
9
 */
 
10
/*
 
11
 * File name:           smyblas2.c
 
12
 * Purpose:
 
13
 *     Level 2 BLAS operations: solves and matvec, written in C.
 
14
 * Note:
 
15
 *     This is only used when the system lacks an efficient BLAS library.
 
16
 */
 
17
 
 
18
/*
 
19
 * Solves a dense UNIT lower triangular system. The unit lower 
 
20
 * triangular matrix is stored in a 2D array M(1:nrow,1:ncol). 
 
21
 * The solution will be returned in the rhs vector.
 
22
 */
 
23
 
 
24
/* local prototypes*/
 
25
void slsolve ( int, int, float *, float *);
 
26
void susolve ( int, int, float *, float *);
 
27
void smatvec ( int, int, int, float *, float *, float *);
 
28
 
 
29
 
 
30
void slsolve ( int ldm, int ncol, float *M, float *rhs )
 
31
{
 
32
    int k;
 
33
    float x0, x1, x2, x3, x4, x5, x6, x7;
 
34
    float *M0;
 
35
    register float *Mki0, *Mki1, *Mki2, *Mki3, *Mki4, *Mki5, *Mki6, *Mki7;
 
36
    register int firstcol = 0;
 
37
 
 
38
    M0 = &M[0];
 
39
 
 
40
    while ( firstcol < ncol - 7 ) { /* Do 8 columns */
 
41
      Mki0 = M0 + 1;
 
42
      Mki1 = Mki0 + ldm + 1;
 
43
      Mki2 = Mki1 + ldm + 1;
 
44
      Mki3 = Mki2 + ldm + 1;
 
45
      Mki4 = Mki3 + ldm + 1;
 
46
      Mki5 = Mki4 + ldm + 1;
 
47
      Mki6 = Mki5 + ldm + 1;
 
48
      Mki7 = Mki6 + ldm + 1;
 
49
 
 
50
      x0 = rhs[firstcol];
 
51
      x1 = rhs[firstcol+1] - x0 * *Mki0++;
 
52
      x2 = rhs[firstcol+2] - x0 * *Mki0++ - x1 * *Mki1++;
 
53
      x3 = rhs[firstcol+3] - x0 * *Mki0++ - x1 * *Mki1++ - x2 * *Mki2++;
 
54
      x4 = rhs[firstcol+4] - x0 * *Mki0++ - x1 * *Mki1++ - x2 * *Mki2++
 
55
                           - x3 * *Mki3++;
 
56
      x5 = rhs[firstcol+5] - x0 * *Mki0++ - x1 * *Mki1++ - x2 * *Mki2++
 
57
                           - x3 * *Mki3++ - x4 * *Mki4++;
 
58
      x6 = rhs[firstcol+6] - x0 * *Mki0++ - x1 * *Mki1++ - x2 * *Mki2++
 
59
                           - x3 * *Mki3++ - x4 * *Mki4++ - x5 * *Mki5++;
 
60
      x7 = rhs[firstcol+7] - x0 * *Mki0++ - x1 * *Mki1++ - x2 * *Mki2++
 
61
                           - x3 * *Mki3++ - x4 * *Mki4++ - x5 * *Mki5++
 
62
                           - x6 * *Mki6++;
 
63
 
 
64
      rhs[++firstcol] = x1;
 
65
      rhs[++firstcol] = x2;
 
66
      rhs[++firstcol] = x3;
 
67
      rhs[++firstcol] = x4;
 
68
      rhs[++firstcol] = x5;
 
69
      rhs[++firstcol] = x6;
 
70
      rhs[++firstcol] = x7;
 
71
      ++firstcol;
 
72
    
 
73
      for (k = firstcol; k < ncol; k++)
 
74
        rhs[k] = rhs[k] - x0 * *Mki0++ - x1 * *Mki1++
 
75
                        - x2 * *Mki2++ - x3 * *Mki3++
 
76
                        - x4 * *Mki4++ - x5 * *Mki5++
 
77
                        - x6 * *Mki6++ - x7 * *Mki7++;
 
78
 
 
79
      M0 += 8 * ldm + 8;
 
80
    }
 
81
 
 
82
    while ( firstcol < ncol - 3 ) { /* Do 4 columns */
 
83
      Mki0 = M0 + 1;
 
84
      Mki1 = Mki0 + ldm + 1;
 
85
      Mki2 = Mki1 + ldm + 1;
 
86
      Mki3 = Mki2 + ldm + 1;
 
87
 
 
88
      x0 = rhs[firstcol];
 
89
      x1 = rhs[firstcol+1] - x0 * *Mki0++;
 
90
      x2 = rhs[firstcol+2] - x0 * *Mki0++ - x1 * *Mki1++;
 
91
      x3 = rhs[firstcol+3] - x0 * *Mki0++ - x1 * *Mki1++ - x2 * *Mki2++;
 
92
 
 
93
      rhs[++firstcol] = x1;
 
94
      rhs[++firstcol] = x2;
 
95
      rhs[++firstcol] = x3;
 
96
      ++firstcol;
 
97
    
 
98
      for (k = firstcol; k < ncol; k++)
 
99
        rhs[k] = rhs[k] - x0 * *Mki0++ - x1 * *Mki1++
 
100
                        - x2 * *Mki2++ - x3 * *Mki3++;
 
101
 
 
102
      M0 += 4 * ldm + 4;
 
103
    }
 
104
 
 
105
    if ( firstcol < ncol - 1 ) { /* Do 2 columns */
 
106
      Mki0 = M0 + 1;
 
107
      Mki1 = Mki0 + ldm + 1;
 
108
 
 
109
      x0 = rhs[firstcol];
 
110
      x1 = rhs[firstcol+1] - x0 * *Mki0++;
 
111
 
 
112
      rhs[++firstcol] = x1;
 
113
      ++firstcol;
 
114
    
 
115
      for (k = firstcol; k < ncol; k++)
 
116
        rhs[k] = rhs[k] - x0 * *Mki0++ - x1 * *Mki1++;
 
117
 
 
118
    }
 
119
    
 
120
}
 
121
 
 
122
/*
 
123
 * Solves a dense upper triangular system. The upper triangular matrix is
 
124
 * stored in a 2-dim array M(1:ldm,1:ncol). The solution will be returned
 
125
 * in the rhs vector.
 
126
 */
 
127
void
 
128
susolve ( ldm, ncol, M, rhs )
 
129
int ldm;        /* in */
 
130
int ncol;       /* in */
 
131
float *M;       /* in */
 
132
float *rhs;     /* modified */
 
133
{
 
134
    float xj;
 
135
    int jcol, j, irow;
 
136
 
 
137
    jcol = ncol - 1;
 
138
 
 
139
    for (j = 0; j < ncol; j++) {
 
140
 
 
141
        xj = rhs[jcol] / M[jcol + jcol*ldm];            /* M(jcol, jcol) */
 
142
        rhs[jcol] = xj;
 
143
        
 
144
        for (irow = 0; irow < jcol; irow++)
 
145
            rhs[irow] -= xj * M[irow + jcol*ldm];       /* M(irow, jcol) */
 
146
 
 
147
        jcol--;
 
148
 
 
149
    }
 
150
}
 
151
 
 
152
 
 
153
/*
 
154
 * Performs a dense matrix-vector multiply: Mxvec = Mxvec + M * vec.
 
155
 * The input matrix is M(1:nrow,1:ncol); The product is returned in Mxvec[].
 
156
 */
 
157
void smatvec ( ldm, nrow, ncol, M, vec, Mxvec )
 
158
 
 
159
int ldm;        /* in -- leading dimension of M */
 
160
int nrow;       /* in */ 
 
161
int ncol;       /* in */
 
162
float *M;       /* in */
 
163
float *vec;     /* in */
 
164
float *Mxvec;   /* in/out */
 
165
 
 
166
{
 
167
    float vi0, vi1, vi2, vi3, vi4, vi5, vi6, vi7;
 
168
    float *M0;
 
169
    register float *Mki0, *Mki1, *Mki2, *Mki3, *Mki4, *Mki5, *Mki6, *Mki7;
 
170
    register int firstcol = 0;
 
171
    int k;
 
172
 
 
173
    M0 = &M[0];
 
174
    while ( firstcol < ncol - 7 ) {     /* Do 8 columns */
 
175
 
 
176
        Mki0 = M0;
 
177
        Mki1 = Mki0 + ldm;
 
178
        Mki2 = Mki1 + ldm;
 
179
        Mki3 = Mki2 + ldm;
 
180
        Mki4 = Mki3 + ldm;
 
181
        Mki5 = Mki4 + ldm;
 
182
        Mki6 = Mki5 + ldm;
 
183
        Mki7 = Mki6 + ldm;
 
184
 
 
185
        vi0 = vec[firstcol++];
 
186
        vi1 = vec[firstcol++];
 
187
        vi2 = vec[firstcol++];
 
188
        vi3 = vec[firstcol++];  
 
189
        vi4 = vec[firstcol++];
 
190
        vi5 = vec[firstcol++];
 
191
        vi6 = vec[firstcol++];
 
192
        vi7 = vec[firstcol++];  
 
193
 
 
194
        for (k = 0; k < nrow; k++) 
 
195
            Mxvec[k] += vi0 * *Mki0++ + vi1 * *Mki1++
 
196
                      + vi2 * *Mki2++ + vi3 * *Mki3++ 
 
197
                      + vi4 * *Mki4++ + vi5 * *Mki5++
 
198
                      + vi6 * *Mki6++ + vi7 * *Mki7++;
 
199
 
 
200
        M0 += 8 * ldm;
 
201
    }
 
202
 
 
203
    while ( firstcol < ncol - 3 ) {     /* Do 4 columns */
 
204
 
 
205
        Mki0 = M0;
 
206
        Mki1 = Mki0 + ldm;
 
207
        Mki2 = Mki1 + ldm;
 
208
        Mki3 = Mki2 + ldm;
 
209
 
 
210
        vi0 = vec[firstcol++];
 
211
        vi1 = vec[firstcol++];
 
212
        vi2 = vec[firstcol++];
 
213
        vi3 = vec[firstcol++];  
 
214
        for (k = 0; k < nrow; k++) 
 
215
            Mxvec[k] += vi0 * *Mki0++ + vi1 * *Mki1++
 
216
                      + vi2 * *Mki2++ + vi3 * *Mki3++ ;
 
217
 
 
218
        M0 += 4 * ldm;
 
219
    }
 
220
 
 
221
    while ( firstcol < ncol ) {         /* Do 1 column */
 
222
 
 
223
        Mki0 = M0;
 
224
        vi0 = vec[firstcol++];
 
225
        for (k = 0; k < nrow; k++)
 
226
            Mxvec[k] += vi0 * *Mki0++;
 
227
 
 
228
        M0 += ldm;
 
229
    }
 
230
        
 
231
}
 
232