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