ergo
template_lapack_lascl.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /* This file belongs to the template_lapack part of the Ergo source
31 * code. The source files in the template_lapack directory are modified
32 * versions of files originally distributed as CLAPACK, see the
33 * Copyright/license notice in the file template_lapack/COPYING.
34 */
35
36
37#ifndef TEMPLATE_LAPACK_LASCL_HEADER
38#define TEMPLATE_LAPACK_LASCL_HEADER
39
40
41template<class Treal>
42int template_lapack_lascl(const char *type__, const integer *kl, const integer *ku,
43 const Treal *cfrom, const Treal *cto, const integer *m, const integer *n,
44 Treal *a, const integer *lda, integer *info)
45{
46/* -- LAPACK auxiliary routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 February 29, 1992
50
51
52 Purpose
53 =======
54
55 DLASCL multiplies the M by N real matrix A by the real scalar
56 CTO/CFROM. This is done without over/underflow as long as the final
57 result CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that
58 A may be full, upper triangular, lower triangular, upper Hessenberg,
59 or banded.
60
61 Arguments
62 =========
63
64 TYPE (input) CHARACTER*1
65 TYPE indices the storage type of the input matrix.
66 = 'G': A is a full matrix.
67 = 'L': A is a lower triangular matrix.
68 = 'U': A is an upper triangular matrix.
69 = 'H': A is an upper Hessenberg matrix.
70 = 'B': A is a symmetric band matrix with lower bandwidth KL
71 and upper bandwidth KU and with the only the lower
72 half stored.
73 = 'Q': A is a symmetric band matrix with lower bandwidth KL
74 and upper bandwidth KU and with the only the upper
75 half stored.
76 = 'Z': A is a band matrix with lower bandwidth KL and upper
77 bandwidth KU.
78
79 KL (input) INTEGER
80 The lower bandwidth of A. Referenced only if TYPE = 'B',
81 'Q' or 'Z'.
82
83 KU (input) INTEGER
84 The upper bandwidth of A. Referenced only if TYPE = 'B',
85 'Q' or 'Z'.
86
87 CFROM (input) DOUBLE PRECISION
88 CTO (input) DOUBLE PRECISION
89 The matrix A is multiplied by CTO/CFROM. A(I,J) is computed
90 without over/underflow if the final result CTO*A(I,J)/CFROM
91 can be represented without over/underflow. CFROM must be
92 nonzero.
93
94 M (input) INTEGER
95 The number of rows of the matrix A. M >= 0.
96
97 N (input) INTEGER
98 The number of columns of the matrix A. N >= 0.
99
100 A (input/output) DOUBLE PRECISION array, dimension (LDA,M)
101 The matrix to be multiplied by CTO/CFROM. See TYPE for the
102 storage type.
103
104 LDA (input) INTEGER
105 The leading dimension of the array A. LDA >= max(1,M).
106
107 INFO (output) INTEGER
108 0 - successful exit
109 <0 - if INFO = -i, the i-th argument had an illegal value.
110
111 =====================================================================
112
113
114 Test the input arguments
115
116 Parameter adjustments */
117 /* System generated locals */
118 integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
119 /* Local variables */
120 logical done;
121 Treal ctoc;
122 integer i__, j;
123 integer itype, k1, k2, k3, k4;
124 Treal cfrom1;
125 Treal cfromc;
126 Treal bignum, smlnum, mul, cto1;
127#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
128
129 a_dim1 = *lda;
130 a_offset = 1 + a_dim1 * 1;
131 a -= a_offset;
132
133 /* Function Body */
134 *info = 0;
135
136 if (template_blas_lsame(type__, "G")) {
137 itype = 0;
138 } else if (template_blas_lsame(type__, "L")) {
139 itype = 1;
140 } else if (template_blas_lsame(type__, "U")) {
141 itype = 2;
142 } else if (template_blas_lsame(type__, "H")) {
143 itype = 3;
144 } else if (template_blas_lsame(type__, "B")) {
145 itype = 4;
146 } else if (template_blas_lsame(type__, "Q")) {
147 itype = 5;
148 } else if (template_blas_lsame(type__, "Z")) {
149 itype = 6;
150 } else {
151 itype = -1;
152 }
153
154 if (itype == -1) {
155 *info = -1;
156 } else if (*cfrom == 0.) {
157 *info = -4;
158 } else if (*m < 0) {
159 *info = -6;
160 } else if (*n < 0 || ( itype == 4 && *n != *m ) || ( itype == 5 && *n != *m ) ) {
161 *info = -7;
162 } else if (itype <= 3 && *lda < maxMACRO(1,*m)) {
163 *info = -9;
164 } else if (itype >= 4) {
165/* Computing MAX */
166 i__1 = *m - 1;
167 if (*kl < 0 || *kl > maxMACRO(i__1,0)) {
168 *info = -2;
169 } else /* if(complicated condition) */ {
170/* Computing MAX */
171 i__1 = *n - 1;
172 if (*ku < 0 || *ku > maxMACRO(i__1,0) || ( (itype == 4 || itype == 5) &&
173 *kl != *ku ) ) {
174 *info = -3;
175 } else if ( ( itype == 4 && *lda < *kl + 1 ) || ( itype == 5 && *lda < *
176 ku + 1 ) || ( itype == 6 && *lda < (*kl << 1) + *ku + 1 ) ) {
177 *info = -9;
178 }
179 }
180 }
181
182 if (*info != 0) {
183 i__1 = -(*info);
184 template_blas_erbla("LASCL ", &i__1);
185 return 0;
186 }
187
188/* Quick return if possible */
189
190 if (*n == 0 || *m == 0) {
191 return 0;
192 }
193
194/* Get machine parameters */
195
196 smlnum = template_lapack_lamch("S", (Treal)0);
197 bignum = 1. / smlnum;
198
199 cfromc = *cfrom;
200 ctoc = *cto;
201
202L10:
203 cfrom1 = cfromc * smlnum;
204 cto1 = ctoc / bignum;
205 if (absMACRO(cfrom1) > absMACRO(ctoc) && ctoc != 0.) {
206 mul = smlnum;
207 done = FALSE_;
208 cfromc = cfrom1;
209 } else if (absMACRO(cto1) > absMACRO(cfromc)) {
210 mul = bignum;
211 done = FALSE_;
212 ctoc = cto1;
213 } else {
214 mul = ctoc / cfromc;
215 done = TRUE_;
216 }
217
218 if (itype == 0) {
219
220/* Full matrix */
221
222 i__1 = *n;
223 for (j = 1; j <= i__1; ++j) {
224 i__2 = *m;
225 for (i__ = 1; i__ <= i__2; ++i__) {
226 a_ref(i__, j) = a_ref(i__, j) * mul;
227/* L20: */
228 }
229/* L30: */
230 }
231
232 } else if (itype == 1) {
233
234/* Lower triangular matrix */
235
236 i__1 = *n;
237 for (j = 1; j <= i__1; ++j) {
238 i__2 = *m;
239 for (i__ = j; i__ <= i__2; ++i__) {
240 a_ref(i__, j) = a_ref(i__, j) * mul;
241/* L40: */
242 }
243/* L50: */
244 }
245
246 } else if (itype == 2) {
247
248/* Upper triangular matrix */
249
250 i__1 = *n;
251 for (j = 1; j <= i__1; ++j) {
252 i__2 = minMACRO(j,*m);
253 for (i__ = 1; i__ <= i__2; ++i__) {
254 a_ref(i__, j) = a_ref(i__, j) * mul;
255/* L60: */
256 }
257/* L70: */
258 }
259
260 } else if (itype == 3) {
261
262/* Upper Hessenberg matrix */
263
264 i__1 = *n;
265 for (j = 1; j <= i__1; ++j) {
266/* Computing MIN */
267 i__3 = j + 1;
268 i__2 = minMACRO(i__3,*m);
269 for (i__ = 1; i__ <= i__2; ++i__) {
270 a_ref(i__, j) = a_ref(i__, j) * mul;
271/* L80: */
272 }
273/* L90: */
274 }
275
276 } else if (itype == 4) {
277
278/* Lower half of a symmetric band matrix */
279
280 k3 = *kl + 1;
281 k4 = *n + 1;
282 i__1 = *n;
283 for (j = 1; j <= i__1; ++j) {
284/* Computing MIN */
285 i__3 = k3, i__4 = k4 - j;
286 i__2 = minMACRO(i__3,i__4);
287 for (i__ = 1; i__ <= i__2; ++i__) {
288 a_ref(i__, j) = a_ref(i__, j) * mul;
289/* L100: */
290 }
291/* L110: */
292 }
293
294 } else if (itype == 5) {
295
296/* Upper half of a symmetric band matrix */
297
298 k1 = *ku + 2;
299 k3 = *ku + 1;
300 i__1 = *n;
301 for (j = 1; j <= i__1; ++j) {
302/* Computing MAX */
303 i__2 = k1 - j;
304 i__3 = k3;
305 for (i__ = maxMACRO(i__2,1); i__ <= i__3; ++i__) {
306 a_ref(i__, j) = a_ref(i__, j) * mul;
307/* L120: */
308 }
309/* L130: */
310 }
311
312 } else if (itype == 6) {
313
314/* Band matrix */
315
316 k1 = *kl + *ku + 2;
317 k2 = *kl + 1;
318 k3 = (*kl << 1) + *ku + 1;
319 k4 = *kl + *ku + 1 + *m;
320 i__1 = *n;
321 for (j = 1; j <= i__1; ++j) {
322/* Computing MAX */
323 i__3 = k1 - j;
324/* Computing MIN */
325 i__4 = k3, i__5 = k4 - j;
326 i__2 = minMACRO(i__4,i__5);
327 for (i__ = maxMACRO(i__3,k2); i__ <= i__2; ++i__) {
328 a_ref(i__, j) = a_ref(i__, j) * mul;
329/* L140: */
330 }
331/* L150: */
332 }
333
334 }
335
336 if (! done) {
337 goto L10;
338 }
339
340 return 0;
341
342/* End of DLASCL */
343
344} /* dlascl_ */
345
346#undef a_ref
347
348
349#endif
int template_blas_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition template_blas_common.cc:46
int integer
Definition template_blas_common.h:40
#define absMACRO(x)
Definition template_blas_common.h:47
#define minMACRO(a, b)
Definition template_blas_common.h:46
#define maxMACRO(a, b)
Definition template_blas_common.h:45
bool logical
Definition template_blas_common.h:41
#define TRUE_
Definition template_lapack_common.h:42
#define FALSE_
Definition template_lapack_common.h:43
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition template_lapack_lamch.h:202
int template_lapack_lascl(const char *type__, const integer *kl, const integer *ku, const Treal *cfrom, const Treal *cto, const integer *m, const integer *n, Treal *a, const integer *lda, integer *info)
Definition template_lapack_lascl.h:42
#define a_ref(a_1, a_2)