1 #include "relapack.h"
2 #if XSYTRF_ALLOW_MALLOC
3 #include <stdlib.h>
4 #endif
5
6 static void RELAPACK_dsytrf_rook_rec(const char *, const blasint *, const blasint *, blasint *,
7 double *, const blasint *, blasint *, double *, const blasint *, blasint *);
8
9
10 /** DSYTRF_ROOK computes the factorization of a real symmetric matrix A using the bounded Bunch-Kaufman ("rook") diagonal pivoting method.
11 *
12 * This routine is functionally equivalent to LAPACK's dsytrf_rook.
13 * For details on its interface, see
14 * http://www.netlib.org/lapack/explore-html/db/df4/dsytrf__rook_8f.html
15 * */
RELAPACK_dsytrf_rook(const char * uplo,const blasint * n,double * A,const blasint * ldA,blasint * ipiv,double * Work,const blasint * lWork,blasint * info)16 void RELAPACK_dsytrf_rook(
17 const char *uplo, const blasint *n,
18 double *A, const blasint *ldA, blasint *ipiv,
19 double *Work, const blasint *lWork, blasint *info
20 ) {
21
22 // Required work size
23 const blasint cleanlWork = *n * (*n / 2);
24 blasint minlWork = cleanlWork;
25 #if XSYTRF_ALLOW_MALLOC
26 minlWork = 1;
27 #endif
28
29 // Check arguments
30 const blasint lower = LAPACK(lsame)(uplo, "L");
31 const blasint upper = LAPACK(lsame)(uplo, "U");
32 *info = 0;
33 if (!lower && !upper)
34 *info = -1;
35 else if (*n < 0)
36 *info = -2;
37 else if (*ldA < MAX(1, *n))
38 *info = -4;
39 else if ((*lWork <1 || *lWork < minlWork) && *lWork != -1)
40 *info = -7;
41 else if (*lWork == -1) {
42 // Work size query
43 *Work = cleanlWork;
44 return;
45 }
46
47 // Ensure Work size
48 double *cleanWork = Work;
49 #if XSYTRF_ALLOW_MALLOC
50 if (!*info && *lWork < cleanlWork) {
51 cleanWork = malloc(cleanlWork * sizeof(double));
52 if (!cleanWork)
53 *info = -7;
54 }
55 #endif
56
57 if (*info) {
58 const blasint minfo = -*info;
59 LAPACK(xerbla)("DSYTRF_ROOK", &minfo, strlen("DSYTRF_ROOK"));
60 return;
61 }
62
63 // Clean char * arguments
64 const char cleanuplo = lower ? 'L' : 'U';
65
66 // Dummy argument
67 blasint nout;
68
69 // Recursive kernel
70 RELAPACK_dsytrf_rook_rec(&cleanuplo, n, n, &nout, A, ldA, ipiv, cleanWork, n, info);
71
72 #if XSYTRF_ALLOW_MALLOC
73 if (cleanWork != Work)
74 free(cleanWork);
75 #endif
76 }
77
78
79 /** dsytrf_rook's recursive compute kernel */
RELAPACK_dsytrf_rook_rec(const char * uplo,const blasint * n_full,const blasint * n,blasint * n_out,double * A,const blasint * ldA,blasint * ipiv,double * Work,const blasint * ldWork,blasint * info)80 static void RELAPACK_dsytrf_rook_rec(
81 const char *uplo, const blasint *n_full, const blasint *n, blasint *n_out,
82 double *A, const blasint *ldA, blasint *ipiv,
83 double *Work, const blasint *ldWork, blasint *info
84 ) {
85
86 // top recursion level?
87 const blasint top = *n_full == *n;
88
89 if (*n <= MAX(CROSSOVER_DSYTRF_ROOK, 3)) {
90 // Unblocked
91 if (top) {
92 LAPACK(dsytf2)(uplo, n, A, ldA, ipiv, info);
93 *n_out = *n;
94 } else
95 RELAPACK_dsytrf_rook_rec2(uplo, n_full, n, n_out, A, ldA, ipiv, Work, ldWork, info);
96 return;
97 }
98
99 blasint info1, info2;
100
101 // Constants
102 const double ONE[] = { 1. };
103 const double MONE[] = { -1. };
104 const blasint iONE[] = { 1 };
105
106 const blasint n_rest = *n_full - *n;
107
108 if (*uplo == 'L') {
109 // Splitting (setup)
110 blasint n1 = DREC_SPLIT(*n);
111 blasint n2 = *n - n1;
112
113 // Work_L *
114 double *const Work_L = Work;
115
116 // recursion(A_L)
117 blasint n1_out;
118 RELAPACK_dsytrf_rook_rec(uplo, n_full, &n1, &n1_out, A, ldA, ipiv, Work_L, ldWork, &info1);
119 n1 = n1_out;
120
121 // Splitting (continued)
122 n2 = *n - n1;
123 const blasint n_full2 = *n_full - n1;
124
125 // * *
126 // A_BL A_BR
127 // A_BL_B A_BR_B
128 double *const A_BL = A + n1;
129 double *const A_BR = A + *ldA * n1 + n1;
130 double *const A_BL_B = A + *n;
131 double *const A_BR_B = A + *ldA * n1 + *n;
132
133 // * *
134 // Work_BL Work_BR
135 // * *
136 // (top recursion level: use Work as Work_BR)
137 double *const Work_BL = Work + n1;
138 double *const Work_BR = top ? Work : Work + *ldWork * n1 + n1;
139 const blasint ldWork_BR = top ? n2 : *ldWork;
140
141 // ipiv_T
142 // ipiv_B
143 blasint *const ipiv_B = ipiv + n1;
144
145 // A_BR = A_BR - A_BL Work_BL'
146 RELAPACK_dgemmt(uplo, "N", "T", &n2, &n1, MONE, A_BL, ldA, Work_BL, ldWork, ONE, A_BR, ldA);
147 BLAS(dgemm)("N", "T", &n_rest, &n2, &n1, MONE, A_BL_B, ldA, Work_BL, ldWork, ONE, A_BR_B, ldA);
148
149 // recursion(A_BR)
150 blasint n2_out;
151 RELAPACK_dsytrf_rook_rec(uplo, &n_full2, &n2, &n2_out, A_BR, ldA, ipiv_B, Work_BR, &ldWork_BR, &info2);
152
153 if (n2_out != n2) {
154 // undo 1 column of updates
155 const blasint n_restp1 = n_rest + 1;
156
157 // last column of A_BR
158 double *const A_BR_r = A_BR + *ldA * n2_out + n2_out;
159
160 // last row of A_BL
161 double *const A_BL_b = A_BL + n2_out;
162
163 // last row of Work_BL
164 double *const Work_BL_b = Work_BL + n2_out;
165
166 // A_BR_r = A_BR_r + A_BL_b Work_BL_b'
167 BLAS(dgemv)("N", &n_restp1, &n1, ONE, A_BL_b, ldA, Work_BL_b, ldWork, ONE, A_BR_r, iONE);
168 }
169 n2 = n2_out;
170
171 // shift pivots
172 blasint i;
173 for (i = 0; i < n2; i++)
174 if (ipiv_B[i] > 0)
175 ipiv_B[i] += n1;
176 else
177 ipiv_B[i] -= n1;
178
179 *info = info1 || info2;
180 *n_out = n1 + n2;
181 } else {
182 // Splitting (setup)
183 blasint n2 = DREC_SPLIT(*n);
184 blasint n1 = *n - n2;
185
186 // * Work_R
187 // (top recursion level: use Work as Work_R)
188 double *const Work_R = top ? Work : Work + *ldWork * n1;
189
190 // recursion(A_R)
191 blasint n2_out;
192 RELAPACK_dsytrf_rook_rec(uplo, n_full, &n2, &n2_out, A, ldA, ipiv, Work_R, ldWork, &info2);
193 const blasint n2_diff = n2 - n2_out;
194 n2 = n2_out;
195
196 // Splitting (continued)
197 n1 = *n - n2;
198 const blasint n_full1 = *n_full - n2;
199
200 // * A_TL_T A_TR_T
201 // * A_TL A_TR
202 // * * *
203 double *const A_TL_T = A + *ldA * n_rest;
204 double *const A_TR_T = A + *ldA * (n_rest + n1);
205 double *const A_TL = A + *ldA * n_rest + n_rest;
206 double *const A_TR = A + *ldA * (n_rest + n1) + n_rest;
207
208 // Work_L *
209 // * Work_TR
210 // * *
211 // (top recursion level: Work_R was Work)
212 double *const Work_L = Work;
213 double *const Work_TR = Work + *ldWork * (top ? n2_diff : n1) + n_rest;
214 const blasint ldWork_L = top ? n1 : *ldWork;
215
216 // A_TL = A_TL - A_TR Work_TR'
217 RELAPACK_dgemmt(uplo, "N", "T", &n1, &n2, MONE, A_TR, ldA, Work_TR, ldWork, ONE, A_TL, ldA);
218 BLAS(dgemm)("N", "T", &n_rest, &n1, &n2, MONE, A_TR_T, ldA, Work_TR, ldWork, ONE, A_TL_T, ldA);
219
220 // recursion(A_TL)
221 blasint n1_out;
222 RELAPACK_dsytrf_rook_rec(uplo, &n_full1, &n1, &n1_out, A, ldA, ipiv, Work_L, &ldWork_L, &info1);
223
224 if (n1_out != n1) {
225 // undo 1 column of updates
226 const blasint n_restp1 = n_rest + 1;
227
228 // A_TL_T_l = A_TL_T_l + A_TR_T Work_TR_t'
229 BLAS(dgemv)("N", &n_restp1, &n2, ONE, A_TR_T, ldA, Work_TR, ldWork, ONE, A_TL_T, iONE);
230 }
231 n1 = n1_out;
232
233 *info = info2 || info1;
234 *n_out = n1 + n2;
235 }
236 }
237