1 /*
2 * @BEGIN LICENSE
3 *
4 * ambit: C++ library for the implementation of tensor product calculations
5 * through a clean, concise user interface.
6 *
7 * Copyright (c) 2014-2017 Ambit developers.
8 *
9 * The copyrights for code used from other parties are included in
10 * the corresponding files.
11 *
12 * This file is part of ambit.
13 *
14 * Ambit is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU Lesser General Public License as published by
16 * the Free Software Foundation, version 3.
17 *
18 * Ambit is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU Lesser General Public License for more details.
22 *
23 * You should have received a copy of the GNU Lesser General Public License along
24 * with ambit; if not, write to the Free Software Foundation, Inc.,
25 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 *
27 * @END LICENSE
28 */
29
30 #include "math.h"
31
32 #include <climits>
33 #include <cmath>
34 #include <stdexcept>
35
36 //#include <FCMangle.h>
37
38 #define FC_GLOBAL(lc, UC) lc ## _
39
40 // => BLAS level 1 <=
41 #define F_DSWAP FC_GLOBAL(dswap, DSWAP)
42 #define F_DAXPY FC_GLOBAL(daxpy, DAXPY)
43 #define F_DCOPY FC_GLOBAL(dcopy, DCOPY)
44 #define F_DROT FC_GLOBAL(drot, DROT)
45 #define F_DSCAL FC_GLOBAL(dscal, DSCAL)
46 #define F_DDOT FC_GLOBAL(ddot, DDOT)
47 #define F_DASUM FC_GLOBAL(dasum, DASUM)
48 #define F_DNRM2 FC_GLOBAL(dnrm2, DNRM2)
49 #define F_IDAMAX FC_GLOBAL(idamax, IDAMAX)
50
51 // => BLAS level 2/3 <=
52 #define F_DGBMV FC_GLOBAL(dgbmv, DGBMV)
53 #define F_DGEMM FC_GLOBAL(dgemm, DGEMM)
54 #define F_DGEMV FC_GLOBAL(dgemv, DGEMV)
55 #define F_DGER FC_GLOBAL(dger, DGER)
56 #define F_DSBMV FC_GLOBAL(dsbmv, DSBMV)
57 #define F_DSPMV FC_GLOBAL(dspmv, DSPMV)
58 #define F_DSPR FC_GLOBAL(dspr, DSPR)
59 #define F_DSPR2 FC_GLOBAL(dspr2, DSPR2)
60 #define F_DSYMM FC_GLOBAL(dsymm, DSYMM)
61 #define F_DSYMV FC_GLOBAL(dsymv, DSYMV)
62 #define F_DSYR FC_GLOBAL(dsyr, DSYR)
63 #define F_DSYR2 FC_GLOBAL(dsyr2, DSYR2)
64 #define F_DSYR2K FC_GLOBAL(dsyr2k, DSYR2K)
65 #define F_DSYRK FC_GLOBAL(dsyrk, DSYRK)
66 #define F_DTBMV FC_GLOBAL(dtbmv, DTBMV)
67 #define F_DTBSV FC_GLOBAL(dtbsv, DTBSV)
68 #define F_DTPMV FC_GLOBAL(dtpmv, DTPMV)
69 #define F_DTPSV FC_GLOBAL(dtpsv, DTPSV)
70 #define F_DTRMM FC_GLOBAL(dtrmm, DTRMM)
71 #define F_DTRMV FC_GLOBAL(dtrmv, DTRMV)
72 #define F_DTRSM FC_GLOBAL(dtrsm, DTRSM)
73 #define F_DTRSV FC_GLOBAL(dtrsv, DTRSV)
74
75 extern "C" {
76
77 // => BLAS level 1 <=
78 extern void F_DSWAP(int *length, double *x, int *incx, double *y, int *inc_y);
79 extern void F_DAXPY(int *length, double *a, double *x, int *inc_x, double *y,
80 int *inc_y);
81 extern void F_DCOPY(int *length, double *x, int *inc_x, double *y, int *inc_y);
82 extern void F_DGEMM(char *transa, char *transb, int *m, int *n, int *k,
83 double *alpha, double *A, int *lda, double *B, int *ldb,
84 double *beta, double *C, int *ldc);
85 extern void F_DSYMM(char *side, char *uplo, int *m, int *n, double *alpha,
86 double *A, int *lda, double *B, int *ldb, double *beta,
87 double *C, int *ldc);
88 extern void F_DROT(int *ntot, double *x, int *incx, double *y, int *incy,
89 double *cotheta, double *sintheta);
90 extern void F_DSCAL(int *n, double *alpha, double *vec, int *inc);
91 extern void F_DGEMV(char *transa, int *m, int *n, double *alpha, double *A,
92 int *lda, double *X, int *inc_x, double *beta, double *Y,
93 int *inc_y);
94 extern void F_DSYMV(char *uplo, int *n, double *alpha, double *A, int *lda,
95 double *X, int *inc_x, double *beta, double *Y, int *inc_y);
96 extern void F_DSPMV(char *uplo, int *n, double *alpha, double *A, double *X,
97 int *inc_x, double *beta, double *Y, int *inc_y);
98 extern double F_DDOT(int *n, double *x, int *incx, double *y, int *incy);
99 extern double F_DNRM2(int *n, double *x, int *incx);
100 extern double F_DASUM(int *n, double *x, int *incx);
101 extern int F_IDAMAX(int *n, double *x, int *incx);
102
103 // => BLAS level 2/3 <=
104 extern void F_DGBMV(char *, int *, int *, int *, int *, double *, double *,
105 int *, double *, int *, double *, double *, int *);
106 extern void F_DGEMM(char *, char *, int *, int *, int *, double *, double *,
107 int *, double *, int *, double *, double *, int *);
108 extern void F_DGEMV(char *, int *, int *, double *, double *, int *, double *,
109 int *, double *, double *, int *);
110 extern void F_DGER(int *, int *, double *, double *, int *, double *, int *,
111 double *, int *);
112 extern void F_DSBMV(char *, int *, int *, double *, double *, int *, double *,
113 int *, double *, double *, int *);
114 extern void F_DSPMV(char *, int *, double *, double *, double *, int *,
115 double *, double *, int *);
116 extern void F_DSPR(char *, int *, double *, double *, int *, double *);
117 extern void F_DSPR2(char *, int *, double *, double *, int *, double *, int *,
118 double *);
119 extern void F_DSYMM(char *, char *, int *, int *, double *, double *, int *,
120 double *, int *, double *, double *, int *);
121 extern void F_DSYMV(char *, int *, double *, double *, int *, double *, int *,
122 double *, double *, int *);
123 extern void F_DSYR(char *, int *, double *, double *, int *, double *, int *);
124 extern void F_DSYR2(char *, int *, double *, double *, int *, double *, int *,
125 double *, int *);
126 extern void F_DSYR2K(char *, char *, int *, int *, double *, double *, int *,
127 double *, int *, double *, double *, int *);
128 extern void F_DSYRK(char *, char *, int *, int *, double *, double *, int *,
129 double *, double *, int *);
130 extern void F_DTBMV(char *, char *, char *, int *, int *, double *, int *,
131 double *, int *);
132 extern void F_DTBSV(char *, char *, char *, int *, int *, double *, int *,
133 double *, int *);
134 extern void F_DTPMV(char *, char *, char *, int *, double *, double *, int *);
135 extern void F_DTPSV(char *, char *, char *, int *, double *, double *, int *);
136 extern void F_DTRMM(char *, char *, char *, char *, int *, int *, double *,
137 double *, int *, double *, int *);
138 extern void F_DTRMV(char *, char *, char *, int *, double *, int *, double *,
139 int *);
140 extern void F_DTRSM(char *, char *, char *, char *, int *, int *, double *,
141 double *, int *, double *, int *);
142 extern void F_DTRSV(char *, char *, char *, int *, double *, int *, double *,
143 int *);
144 }
145
146 namespace ambit
147 {
148
149 /**
150 * Swaps a vector with another vector.
151 *
152 * @param length Specifies the number of elements in vectors x and y.
153 * @param x Array, DIMENSION at least (1 + (n-1)*abs(incx)).
154 * @param inc_x Specifies the increment for the elements of x.
155 * @param y Array, DIMENSION at least (1 + (n-1)*abs(incy)).
156 * @param inc_y Specifies the increment for the elements of y.
157 *
158 */
C_DSWAP(unsigned long int length,double * x,int inc_x,double * y,int inc_y)159 void C_DSWAP(unsigned long int length, double *x, int inc_x, double *y,
160 int inc_y)
161 {
162 int big_blocks = (int)(length / INT_MAX);
163 int small_size = (int)(length % INT_MAX);
164 for (int block = 0; block <= big_blocks; block++)
165 {
166 double *x_s = &x[block * inc_x * (unsigned long int)INT_MAX];
167 double *y_s = &y[block * inc_y * (unsigned long int)INT_MAX];
168 signed int length_s = (block == big_blocks) ? small_size : INT_MAX;
169 ::F_DSWAP(&length_s, x_s, &inc_x, y_s, &inc_y);
170 }
171 }
172
173 /*!
174 * This function performs y = a * x + y.
175 *
176 * Steps every inc_x in x and every inc_y in y (normally both 1).
177 *
178 * \param length length of arrays
179 * \param a scalar a to multiply vector x
180 * \param x vector x
181 * \param inc_x how many places to skip to get to next element in x
182 * \param y vector y
183 * \param inc_y how many places to skip to get to next element in y
184 *
185 */
C_DAXPY(unsigned long int length,double a,double * x,int inc_x,double * y,int inc_y)186 void C_DAXPY(unsigned long int length, double a, double *x, int inc_x,
187 double *y, int inc_y)
188 {
189 int big_blocks = (int)(length / INT_MAX);
190 int small_size = (int)(length % INT_MAX);
191 for (int block = 0; block <= big_blocks; block++)
192 {
193 double *x_s = &x[block * inc_x * (unsigned long int)INT_MAX];
194 double *y_s = &y[block * inc_y * (unsigned long int)INT_MAX];
195 signed int length_s = (block == big_blocks) ? small_size : INT_MAX;
196 ::F_DAXPY(&length_s, &a, x_s, &inc_x, y_s, &inc_y);
197 }
198 }
199
200 /*!
201 * This function copies x into y.
202 *
203 * Steps every inc_x in x and every inc_y in y (normally both 1).
204 *
205 * \param length = length of array
206 * \param x = vector x
207 * \param inc_x = how many places to skip to get to next element in x
208 * \param y = vector y
209 * \param inc_y = how many places to skip to get to next element in y
210 *
211 */
C_DCOPY(unsigned long int length,double * x,int inc_x,double * y,int inc_y)212 void C_DCOPY(unsigned long int length, double *x, int inc_x, double *y,
213 int inc_y)
214 {
215 int big_blocks = (int)(length / INT_MAX);
216 int small_size = (int)(length % INT_MAX);
217 for (int block = 0; block <= big_blocks; block++)
218 {
219 double *x_s = &x[block * inc_x * (unsigned long int)INT_MAX];
220 double *y_s = &y[block * inc_y * (unsigned long int)INT_MAX];
221 signed int length_s = (block == big_blocks) ? small_size : INT_MAX;
222 ::F_DCOPY(&length_s, x_s, &inc_x, y_s, &inc_y);
223 }
224 }
225
226 /*!
227 * This function scales a vector by a real scalar.
228 *
229 * \param length length of array
230 * \param alpha scale factor
231 * \param vec vector to scale
232 * \param inc how many places to skip to get to next element in vec
233 *
234 */
C_DSCAL(unsigned long int length,double alpha,double * vec,int inc)235 void C_DSCAL(unsigned long int length, double alpha, double *vec, int inc)
236 {
237 int big_blocks = (int)(length / INT_MAX);
238 int small_size = (int)(length % INT_MAX);
239 for (int block = 0; block <= big_blocks; block++)
240 {
241 double *vec_s = &vec[block * inc * (unsigned long int)INT_MAX];
242 signed int length_s = (block == big_blocks) ? small_size : INT_MAX;
243 ::F_DSCAL(&length_s, &alpha, vec_s, &inc);
244 }
245 }
246
247 /*!
248 *Calculates a plane Givens rotation for vectors x, y and
249 * angle theta. x = x*cos + y*sin, y = -x*sin + y*cos.
250 *
251 * \param x vector x
252 * \param y vector Y
253 * \param length length of x,y
254 * \param inc_x how many places to skip to get to the next element of x
255 * \param inc_y how many places to skip to get to the next element of y
256 *
257 */
C_DROT(unsigned long int length,double * x,int inc_x,double * y,int inc_y,double costheta,double sintheta)258 void C_DROT(unsigned long int length, double *x, int inc_x, double *y,
259 int inc_y, double costheta, double sintheta)
260 {
261 int big_blocks = (int)(length / INT_MAX);
262 int small_size = (int)(length % INT_MAX);
263 for (int block = 0; block <= big_blocks; block++)
264 {
265 double *x_s = &x[block * inc_x * (unsigned long int)INT_MAX];
266 double *y_s = &y[block * inc_y * (unsigned long int)INT_MAX];
267 signed int length_s = (block == big_blocks) ? small_size : INT_MAX;
268 ::F_DROT(&length_s, x_s, &inc_x, y_s, &inc_y, &costheta, &sintheta);
269 }
270 }
271
272 /*!
273 * This function returns the dot product of two vectors, x and y.
274 *
275 * \param length Number of elements in x and y.
276 * \param x A pointer to the beginning of the data in x.
277 * Must be of at least length (1+(N-1)*abs(inc_x).
278 * \param inc_x how many places to skip to get to next element in x
279 * \param y A pointer to the beginning of the data in y.
280 * \param inc_y how many places to skip to get to next element in y
281 *
282 * @return the dot product
283 *
284 */
285
C_DDOT(unsigned long int length,double * x,int inc_x,double * y,int inc_y)286 double C_DDOT(unsigned long int length, double *x, int inc_x, double *y,
287 int inc_y)
288 {
289 if (length == 0)
290 return 0.0;
291
292 double reg = 0.0;
293
294 int big_blocks = (int)(length / INT_MAX);
295 int small_size = (int)(length % INT_MAX);
296 for (int block = 0; block <= big_blocks; block++)
297 {
298 double *x_s = &x[block * inc_x * (unsigned long int)INT_MAX];
299 double *y_s = &y[block * inc_y * (unsigned long int)INT_MAX];
300 signed int length_s = (block == big_blocks) ? small_size : INT_MAX;
301 reg += ::F_DDOT(&length_s, x_s, &inc_x, y_s, &inc_y);
302 }
303
304 return reg;
305 }
306 /*!
307 * This function returns the square of the norm of this vector.
308 *
309 * \param length Number of elements in x.
310 * \param x A pointer to the beginning of the data in x.
311 * Must be of at least length (1+(N-1)*abs(inc_x).
312 * \param inc_x how many places to skip to get to next element in x
313 *
314 * @return the norm squared product
315 *
316 */
317
C_DNRM2(unsigned long int length,double * x,int inc_x)318 double C_DNRM2(unsigned long int length, double *x, int inc_x)
319 {
320 if (length == 0)
321 return 0.0;
322
323 double reg = 0.0;
324
325 int big_blocks = (int)(length / INT_MAX);
326 int small_size = (int)(length % INT_MAX);
327 for (int block = 0; block <= big_blocks; block++)
328 {
329 double *x_s = &x[block * inc_x * (unsigned long int)INT_MAX];
330 signed int length_s = (block == big_blocks) ? small_size : INT_MAX;
331 reg += ::F_DNRM2(&length_s, x_s, &inc_x);
332 }
333
334 return reg;
335 }
336 /*!
337 * This function returns the sum of the absolute value of this vector.
338 *
339 * \param length Number of elements in x.
340 * \param x A pointer to the beginning of the data in x.
341 * Must be of at least length (1+(N-1)*abs(inc_x).
342 * \param inc_x how many places to skip to get to next element in x
343 *
344 * @return the sum of the absolute value
345 *
346 */
347
C_DASUM(unsigned long int length,double * x,int inc_x)348 double C_DASUM(unsigned long int length, double *x, int inc_x)
349 {
350 if (length == 0)
351 return 0.0;
352
353 double reg = 0.0;
354
355 int big_blocks = (int)(length / INT_MAX);
356 int small_size = (int)(length % INT_MAX);
357 for (int block = 0; block <= big_blocks; block++)
358 {
359 double *x_s = &x[block * inc_x * (unsigned long int)INT_MAX];
360 signed int length_s = (block == big_blocks) ? small_size : INT_MAX;
361 reg += ::F_DASUM(&length_s, x_s, &inc_x);
362 }
363
364 return reg;
365 }
366
367 /*!
368 * This function returns the index of the largest absolute value compoment of
369 * this vector.
370 *
371 * \param length Number of elements in x.
372 * \param x A pointer to the beginning of the data in x.
373 * Must be of at least length (1+(N-1)*abs(inc_x).
374 * \param inc_x how many places to skip to get to next element in x
375 *
376 * @return the index of the largest absolute value
377 *
378 */
379
C_IDAMAX(unsigned long int length,double * x,int inc_x)380 unsigned long int C_IDAMAX(unsigned long int length, double *x, int inc_x)
381 {
382 if (length == 0)
383 return 0L;
384
385 unsigned long int reg = 0L;
386 unsigned long int reg2 = 0L;
387
388 int big_blocks = (int)(length / INT_MAX);
389 int small_size = (int)(length % INT_MAX);
390 for (int block = 0; block <= big_blocks; block++)
391 {
392 double *x_s = &x[block * inc_x * (unsigned long int)INT_MAX];
393 signed int length_s = (block == big_blocks) ? small_size : INT_MAX;
394 reg2 = ::F_IDAMAX(&length_s, x_s, &inc_x) +
395 block * inc_x * (unsigned long int)INT_MAX;
396 if (fabs(x[reg]) > fabs(x[reg2]))
397 reg = reg2;
398 }
399
400 return reg;
401 }
402
403 /**
404 * Purpose
405 * =======
406 *
407 * DGBMV performs one of the matrix-vector operations
408 *
409 * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
410 *
411 * where alpha and beta are scalars, x and y are vectors and A is an
412 * m by n band matrix, with kl sub-diagonals and ku super-diagonals.
413 *
414 * Arguments
415 * ==========
416 *
417 * TRANS - CHARACTER*1.
418 * On entry, TRANS specifies the operation to be performed as
419 * follows:
420 *
421 * TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
422 *
423 * TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
424 *
425 * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
426 *
427 * Unchanged on exit.
428 *
429 * M - INTEGER.
430 * On entry, M specifies the number of rows of the matrix A.
431 * M must be at least zero.
432 * Unchanged on exit.
433 *
434 * N - INTEGER.
435 * On entry, N specifies the number of columns of the matrix A.
436 * N must be at least zero.
437 * Unchanged on exit.
438 *
439 * KL - INTEGER.
440 * On entry, KL specifies the number of sub-diagonals of the
441 * matrix A. KL must satisfy 0 .le. KL.
442 * Unchanged on exit.
443 *
444 * KU - INTEGER.
445 * On entry, KU specifies the number of super-diagonals of the
446 * matrix A. KU must satisfy 0 .le. KU.
447 * Unchanged on exit.
448 *
449 * ALPHA - DOUBLE PRECISION.
450 * On entry, ALPHA specifies the scalar alpha.
451 * Unchanged on exit.
452 *
453 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
454 * Before entry, the leading ( kl + ku + 1 ) by n part of the
455 * array A must contain the matrix of coefficients, supplied
456 * column by column, with the leading diagonal of the matrix in
457 * row ( ku + 1 ) of the array, the first super-diagonal
458 * starting at position 2 in row ku, the first sub-diagonal
459 * starting at position 1 in row ( ku + 2 ), and so on.
460 * Elements in the array A that do not correspond to elements
461 * in the band matrix (such as the top left ku by ku triangle)
462 * are not referenced.
463 * The following program segment will transfer a band matrix
464 * from conventional full matrix storage to band storage:
465 *
466 * DO 20, J = 1, N
467 * K = KU + 1 - J
468 * DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
469 * A( K + I, J ) = matrix( I, J )
470 * 10 CONTINUE
471 * 20 CONTINUE
472 *
473 * Unchanged on exit.
474 *
475 * LDA - INTEGER.
476 * On entry, LDA specifies the first dimension of A as declared
477 * in the calling (sub) program. LDA must be at least
478 * ( kl + ku + 1 ).
479 * Unchanged on exit.
480 *
481 * X - DOUBLE PRECISION array of DIMENSION at least
482 * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
483 * and at least
484 * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
485 * Before entry, the incremented array X must contain the
486 * vector x.
487 * Unchanged on exit.
488 *
489 * INCX - INTEGER.
490 * On entry, INCX specifies the increment for the elements of
491 * X. INCX must not be zero.
492 * Unchanged on exit.
493 *
494 * BETA - DOUBLE PRECISION.
495 * On entry, BETA specifies the scalar beta. When BETA is
496 * supplied as zero then Y need not be set on input.
497 * Unchanged on exit.
498 *
499 * Y - DOUBLE PRECISION array of DIMENSION at least
500 * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
501 * and at least
502 * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
503 * Before entry, the incremented array Y must contain the
504 * vector y. On exit, Y is overwritten by the updated vector y.
505 *
506 * INCY - INTEGER.
507 * On entry, INCY specifies the increment for the elements of
508 * Y. INCY must not be zero.
509 * Unchanged on exit.
510 *
511 *
512 * Level 2 Blas routine.
513 *
514 * -- Written on 22-October-1986.
515 * Jack Dongarra, Argonne National Lab.
516 * Jeremy Du Croz, Nag Central Office.
517 * Sven Hammarling, Nag Central Office.
518 * Richard Hanson, Sandia National Labs.
519 *
520 * .. Parameters ..
521 **/
C_DGBMV(char trans,int m,int n,int kl,int ku,double alpha,double * a,int lda,double * x,int incx,double beta,double * y,int incy)522 void C_DGBMV(char trans, int m, int n, int kl, int ku, double alpha, double *a,
523 int lda, double *x, int incx, double beta, double *y, int incy)
524 {
525 if (m == 0 || n == 0)
526 return;
527 if (trans == 'N' || trans == 'n')
528 trans = 'T';
529 else if (trans == 'T' || trans == 't')
530 trans = 'N';
531 else
532 throw std::invalid_argument("C_DGBMV trans argument is invalid.");
533 ::F_DGBMV(&trans, &n, &m, &ku, &kl, &alpha, a, &lda, x, &incx, &beta, y,
534 &incy);
535 }
536
537 /**
538 * Purpose
539 * =======
540 *
541 * DGEMM performs one of the matrix-matrix operations
542 *
543 * C := alpha*op( A )*op( B ) + beta*C,
544 *
545 * where op( X ) is one of
546 *
547 * op( X ) = X or op( X ) = X',
548 *
549 * alpha and beta are scalars, and A, B and C are matrices, with op( A )
550 * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
551 *
552 * Arguments
553 * ==========
554 *
555 * TRANSA - CHARACTER*1.
556 * On entry, TRANSA specifies the form of op( A ) to be used in
557 * the matrix multiplication as follows:
558 *
559 * TRANSA = 'N' or 'n', op( A ) = A.
560 *
561 * TRANSA = 'T' or 't', op( A ) = A'.
562 *
563 * TRANSA = 'C' or 'c', op( A ) = A'.
564 *
565 * Unchanged on exit.
566 *
567 * TRANSB - CHARACTER*1.
568 * On entry, TRANSB specifies the form of op( B ) to be used in
569 * the matrix multiplication as follows:
570 *
571 * TRANSB = 'N' or 'n', op( B ) = B.
572 *
573 * TRANSB = 'T' or 't', op( B ) = B'.
574 *
575 * TRANSB = 'C' or 'c', op( B ) = B'.
576 *
577 * Unchanged on exit.
578 *
579 * M - INTEGER.
580 * On entry, M specifies the number of rows of the matrix
581 * op( A ) and of the matrix C. M must be at least zero.
582 * Unchanged on exit.
583 *
584 * N - INTEGER.
585 * On entry, N specifies the number of columns of the matrix
586 * op( B ) and the number of columns of the matrix C. N must be
587 * at least zero.
588 * Unchanged on exit.
589 *
590 * K - INTEGER.
591 * On entry, K specifies the number of columns of the matrix
592 * op( A ) and the number of rows of the matrix op( B ). K must
593 * be at least zero.
594 * Unchanged on exit.
595 *
596 * ALPHA - DOUBLE PRECISION.
597 * On entry, ALPHA specifies the scalar alpha.
598 * Unchanged on exit.
599 *
600 * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
601 * k when TRANSA = 'N' or 'n', and is m otherwise.
602 * Before entry with TRANSA = 'N' or 'n', the leading m by k
603 * part of the array A must contain the matrix A, otherwise
604 * the leading k by m part of the array A must contain the
605 * matrix A.
606 * Unchanged on exit.
607 *
608 * LDA - INTEGER.
609 * On entry, LDA specifies the first dimension of A as declared
610 * in the calling (sub) program. When TRANSA = 'N' or 'n' then
611 * LDA must be at least max( 1, m ), otherwise LDA must be at
612 * least max( 1, k ).
613 * Unchanged on exit.
614 *
615 * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
616 * n when TRANSB = 'N' or 'n', and is k otherwise.
617 * Before entry with TRANSB = 'N' or 'n', the leading k by n
618 * part of the array B must contain the matrix B, otherwise
619 * the leading n by k part of the array B must contain the
620 * matrix B.
621 * Unchanged on exit.
622 *
623 * LDB - INTEGER.
624 * On entry, LDB specifies the first dimension of B as declared
625 * in the calling (sub) program. When TRANSB = 'N' or 'n' then
626 * LDB must be at least max( 1, k ), otherwise LDB must be at
627 * least max( 1, n ).
628 * Unchanged on exit.
629 *
630 * BETA - DOUBLE PRECISION.
631 * On entry, BETA specifies the scalar beta. When BETA is
632 * supplied as zero then C need not be set on input.
633 * Unchanged on exit.
634 *
635 * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
636 * Before entry, the leading m by n part of the array C must
637 * contain the matrix C, except when beta is zero, in which
638 * case C need not be set on entry.
639 * On exit, the array C is overwritten by the m by n matrix
640 * ( alpha*op( A )*op( B ) + beta*C ).
641 *
642 * LDC - INTEGER.
643 * On entry, LDC specifies the first dimension of C as declared
644 * in the calling (sub) program. LDC must be at least
645 * max( 1, m ).
646 * Unchanged on exit.
647 *
648 *
649 * Level 3 Blas routine.
650 *
651 * -- Written on 8-February-1989.
652 * Jack Dongarra, Argonne National Laboratory.
653 * Iain Duff, AERE Harwell.
654 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
655 * Sven Hammarling, Numerical Algorithms Group Ltd.
656 *
657 *
658 * .. External Functions ..
659 **/
C_DGEMM(char transa,char transb,int m,int n,int k,double alpha,double * a,int lda,double * b,int ldb,double beta,double * c,int ldc)660 void C_DGEMM(char transa, char transb, int m, int n, int k, double alpha,
661 double *a, int lda, double *b, int ldb, double beta, double *c,
662 int ldc)
663 {
664 if (m == 0 || n == 0 || k == 0)
665 return;
666 ::F_DGEMM(&transb, &transa, &n, &m, &k, &alpha, b, &ldb, a, &lda, &beta, c,
667 &ldc);
668 }
669
670 /**
671 * Purpose
672 * =======
673 *
674 * DGEMV performs one of the matrix-vector operations
675 *
676 * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
677 *
678 * where alpha and beta are scalars, x and y are vectors and A is an
679 * m by n matrix.
680 *
681 * Arguments
682 * ==========
683 *
684 * TRANS - CHARACTER*1.
685 * On entry, TRANS specifies the operation to be performed as
686 * follows:
687 *
688 * TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
689 *
690 * TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
691 *
692 * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
693 *
694 * Unchanged on exit.
695 *
696 * M - INTEGER.
697 * On entry, M specifies the number of rows of the matrix A.
698 * M must be at least zero.
699 * Unchanged on exit.
700 *
701 * N - INTEGER.
702 * On entry, N specifies the number of columns of the matrix A.
703 * N must be at least zero.
704 * Unchanged on exit.
705 *
706 * ALPHA - DOUBLE PRECISION.
707 * On entry, ALPHA specifies the scalar alpha.
708 * Unchanged on exit.
709 *
710 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
711 * Before entry, the leading m by n part of the array A must
712 * contain the matrix of coefficients.
713 * Unchanged on exit.
714 *
715 * LDA - INTEGER.
716 * On entry, LDA specifies the first dimension of A as declared
717 * in the calling (sub) program. LDA must be at least
718 * max( 1, m ).
719 * Unchanged on exit.
720 *
721 * X - DOUBLE PRECISION array of DIMENSION at least
722 * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
723 * and at least
724 * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
725 * Before entry, the incremented array X must contain the
726 * vector x.
727 * Unchanged on exit.
728 *
729 * INCX - INTEGER.
730 * On entry, INCX specifies the increment for the elements of
731 * X. INCX must not be zero.
732 * Unchanged on exit.
733 *
734 * BETA - DOUBLE PRECISION.
735 * On entry, BETA specifies the scalar beta. When BETA is
736 * supplied as zero then Y need not be set on input.
737 * Unchanged on exit.
738 *
739 * Y - DOUBLE PRECISION array of DIMENSION at least
740 * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
741 * and at least
742 * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
743 * Before entry with BETA non-zero, the incremented array Y
744 * must contain the vector y. On exit, Y is overwritten by the
745 * updated vector y.
746 *
747 * INCY - INTEGER.
748 * On entry, INCY specifies the increment for the elements of
749 * Y. INCY must not be zero.
750 * Unchanged on exit.
751 *
752 *
753 * Level 2 Blas routine.
754 *
755 * -- Written on 22-October-1986.
756 * Jack Dongarra, Argonne National Lab.
757 * Jeremy Du Croz, Nag Central Office.
758 * Sven Hammarling, Nag Central Office.
759 * Richard Hanson, Sandia National Labs.
760 *
761 *
762 * .. Parameters ..
763 **/
C_DGEMV(char trans,int m,int n,double alpha,double * a,int lda,double * x,int incx,double beta,double * y,int incy)764 void C_DGEMV(char trans, int m, int n, double alpha, double *a, int lda,
765 double *x, int incx, double beta, double *y, int incy)
766 {
767 if (m == 0 || n == 0)
768 return;
769 if (trans == 'N' || trans == 'n')
770 trans = 'T';
771 else if (trans == 'T' || trans == 't')
772 trans = 'N';
773 else
774 throw std::invalid_argument("C_DGEMV trans argument is invalid.");
775 ::F_DGEMV(&trans, &n, &m, &alpha, a, &lda, x, &incx, &beta, y, &incy);
776 }
777
778 /**
779 * Purpose
780 * =======
781 *
782 * DGER performs the rank 1 operation
783 *
784 * A := alpha*x*y' + A,
785 *
786 * where alpha is a scalar, x is an m element vector, y is an n element
787 * vector and A is an m by n matrix.
788 *
789 * Arguments
790 * ==========
791 *
792 * M - INTEGER.
793 * On entry, M specifies the number of rows of the matrix A.
794 * M must be at least zero.
795 * Unchanged on exit.
796 *
797 * N - INTEGER.
798 * On entry, N specifies the number of columns of the matrix A.
799 * N must be at least zero.
800 * Unchanged on exit.
801 *
802 * ALPHA - DOUBLE PRECISION.
803 * On entry, ALPHA specifies the scalar alpha.
804 * Unchanged on exit.
805 *
806 * X - DOUBLE PRECISION array of dimension at least
807 * ( 1 + ( m - 1 )*abs( INCX ) ).
808 * Before entry, the incremented array X must contain the m
809 * element vector x.
810 * Unchanged on exit.
811 *
812 * INCX - INTEGER.
813 * On entry, INCX specifies the increment for the elements of
814 * X. INCX must not be zero.
815 * Unchanged on exit.
816 *
817 * Y - DOUBLE PRECISION array of dimension at least
818 * ( 1 + ( n - 1 )*abs( INCY ) ).
819 * Before entry, the incremented array Y must contain the n
820 * element vector y.
821 * Unchanged on exit.
822 *
823 * INCY - INTEGER.
824 * On entry, INCY specifies the increment for the elements of
825 * Y. INCY must not be zero.
826 * Unchanged on exit.
827 *
828 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
829 * Before entry, the leading m by n part of the array A must
830 * contain the matrix of coefficients. On exit, A is
831 * overwritten by the updated matrix.
832 *
833 * LDA - INTEGER.
834 * On entry, LDA specifies the first dimension of A as declared
835 * in the calling (sub) program. LDA must be at least
836 * max( 1, m ).
837 * Unchanged on exit.
838 *
839 *
840 * Level 2 Blas routine.
841 *
842 * -- Written on 22-October-1986.
843 * Jack Dongarra, Argonne National Lab.
844 * Jeremy Du Croz, Nag Central Office.
845 * Sven Hammarling, Nag Central Office.
846 * Richard Hanson, Sandia National Labs.
847 *
848 *
849 * .. Parameters ..
850 **/
C_DGER(int m,int n,double alpha,double * x,int incx,double * y,int incy,double * a,int lda)851 void C_DGER(int m, int n, double alpha, double *x, int incx, double *y,
852 int incy, double *a, int lda)
853 {
854 if (m == 0 || n == 0)
855 return;
856 ::F_DGER(&n, &m, &alpha, y, &incy, x, &incx, a, &lda);
857 }
858
859 /**
860 * Purpose
861 * =======
862 *
863 * DSBMV performs the matrix-vector operation
864 *
865 * y := alpha*A*x + beta*y,
866 *
867 * where alpha and beta are scalars, x and y are n element vectors and
868 * A is an n by n symmetric band matrix, with k super-diagonals.
869 *
870 * Arguments
871 * ==========
872 *
873 * UPLO - CHARACTER*1.
874 * On entry, UPLO specifies whether the upper or lower
875 * triangular part of the band matrix A is being supplied as
876 * follows:
877 *
878 * UPLO = 'U' or 'u' The upper triangular part of A is
879 * being supplied.
880 *
881 * UPLO = 'L' or 'l' The lower triangular part of A is
882 * being supplied.
883 *
884 * Unchanged on exit.
885 *
886 * N - INTEGER.
887 * On entry, N specifies the order of the matrix A.
888 * N must be at least zero.
889 * Unchanged on exit.
890 *
891 * K - INTEGER.
892 * On entry, K specifies the number of super-diagonals of the
893 * matrix A. K must satisfy 0 .le. K.
894 * Unchanged on exit.
895 *
896 * ALPHA - DOUBLE PRECISION.
897 * On entry, ALPHA specifies the scalar alpha.
898 * Unchanged on exit.
899 *
900 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
901 * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
902 * by n part of the array A must contain the upper triangular
903 * band part of the symmetric matrix, supplied column by
904 * column, with the leading diagonal of the matrix in row
905 * ( k + 1 ) of the array, the first super-diagonal starting at
906 * position 2 in row k, and so on. The top left k by k triangle
907 * of the array A is not referenced.
908 * The following program segment will transfer the upper
909 * triangular part of a symmetric band matrix from conventional
910 * full matrix storage to band storage:
911 *
912 * DO 20, J = 1, N
913 * M = K + 1 - J
914 * DO 10, I = MAX( 1, J - K ), J
915 * A( M + I, J ) = matrix( I, J )
916 * 10 CONTINUE
917 * 20 CONTINUE
918 *
919 * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
920 * by n part of the array A must contain the lower triangular
921 * band part of the symmetric matrix, supplied column by
922 * column, with the leading diagonal of the matrix in row 1 of
923 * the array, the first sub-diagonal starting at position 1 in
924 * row 2, and so on. The bottom right k by k triangle of the
925 * array A is not referenced.
926 * The following program segment will transfer the lower
927 * triangular part of a symmetric band matrix from conventional
928 * full matrix storage to band storage:
929 *
930 * DO 20, J = 1, N
931 * M = 1 - J
932 * DO 10, I = J, MIN( N, J + K )
933 * A( M + I, J ) = matrix( I, J )
934 * 10 CONTINUE
935 * 20 CONTINUE
936 *
937 * Unchanged on exit.
938 *
939 * LDA - INTEGER.
940 * On entry, LDA specifies the first dimension of A as declared
941 * in the calling (sub) program. LDA must be at least
942 * ( k + 1 ).
943 * Unchanged on exit.
944 *
945 * X - DOUBLE PRECISION array of DIMENSION at least
946 * ( 1 + ( n - 1 )*abs( INCX ) ).
947 * Before entry, the incremented array X must contain the
948 * vector x.
949 * Unchanged on exit.
950 *
951 * INCX - INTEGER.
952 * On entry, INCX specifies the increment for the elements of
953 * X. INCX must not be zero.
954 * Unchanged on exit.
955 *
956 * BETA - DOUBLE PRECISION.
957 * On entry, BETA specifies the scalar beta.
958 * Unchanged on exit.
959 *
960 * Y - DOUBLE PRECISION array of DIMENSION at least
961 * ( 1 + ( n - 1 )*abs( INCY ) ).
962 * Before entry, the incremented array Y must contain the
963 * vector y. On exit, Y is overwritten by the updated vector y.
964 *
965 * INCY - INTEGER.
966 * On entry, INCY specifies the increment for the elements of
967 * Y. INCY must not be zero.
968 * Unchanged on exit.
969 *
970 *
971 * Level 2 Blas routine.
972 *
973 * -- Written on 22-October-1986.
974 * Jack Dongarra, Argonne National Lab.
975 * Jeremy Du Croz, Nag Central Office.
976 * Sven Hammarling, Nag Central Office.
977 * Richard Hanson, Sandia National Labs.
978 *
979 *
980 * .. Parameters ..
981 **/
C_DSBMV(char uplo,int n,int k,double alpha,double * a,int lda,double * x,int incx,double beta,double * y,int incy)982 void C_DSBMV(char uplo, int n, int k, double alpha, double *a, int lda,
983 double *x, int incx, double beta, double *y, int incy)
984 {
985 if (n == 0)
986 return;
987 if (uplo == 'U' || uplo == 'u')
988 uplo = 'L';
989 else if (uplo == 'L' || uplo == 'l')
990 uplo = 'U';
991 else
992 throw std::invalid_argument("C_DSBMV uplo argument is invalid.");
993 ::F_DSBMV(&uplo, &n, &k, &alpha, a, &lda, x, &incx, &beta, y, &incy);
994 }
995
996 /**
997 * Purpose
998 * =======
999 *
1000 * DSPMV performs the matrix-vector operation
1001 *
1002 * y := alpha*A*x + beta*y,
1003 *
1004 * where alpha and beta are scalars, x and y are n element vectors and
1005 * A is an n by n symmetric matrix, supplied in packed form.
1006 *
1007 * Arguments
1008 * ==========
1009 *
1010 * UPLO - CHARACTER*1.
1011 * On entry, UPLO specifies whether the upper or lower
1012 * triangular part of the matrix A is supplied in the packed
1013 * array AP as follows:
1014 *
1015 * UPLO = 'U' or 'u' The upper triangular part of A is
1016 * supplied in AP.
1017 *
1018 * UPLO = 'L' or 'l' The lower triangular part of A is
1019 * supplied in AP.
1020 *
1021 * Unchanged on exit.
1022 *
1023 * N - INTEGER.
1024 * On entry, N specifies the order of the matrix A.
1025 * N must be at least zero.
1026 * Unchanged on exit.
1027 *
1028 * ALPHA - DOUBLE PRECISION.
1029 * On entry, ALPHA specifies the scalar alpha.
1030 * Unchanged on exit.
1031 *
1032 * AP - DOUBLE PRECISION array of DIMENSION at least
1033 * ( ( n*( n + 1 ) )/2 ).
1034 * Before entry with UPLO = 'U' or 'u', the array AP must
1035 * contain the upper triangular part of the symmetric matrix
1036 * packed sequentially, column by column, so that AP( 1 )
1037 * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
1038 * and a( 2, 2 ) respectively, and so on.
1039 * Before entry with UPLO = 'L' or 'l', the array AP must
1040 * contain the lower triangular part of the symmetric matrix
1041 * packed sequentially, column by column, so that AP( 1 )
1042 * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
1043 * and a( 3, 1 ) respectively, and so on.
1044 * Unchanged on exit.
1045 *
1046 * X - DOUBLE PRECISION array of dimension at least
1047 * ( 1 + ( n - 1 )*abs( INCX ) ).
1048 * Before entry, the incremented array X must contain the n
1049 * element vector x.
1050 * Unchanged on exit.
1051 *
1052 * INCX - INTEGER.
1053 * On entry, INCX specifies the increment for the elements of
1054 * X. INCX must not be zero.
1055 * Unchanged on exit.
1056 *
1057 * BETA - DOUBLE PRECISION.
1058 * On entry, BETA specifies the scalar beta. When BETA is
1059 * supplied as zero then Y need not be set on input.
1060 * Unchanged on exit.
1061 *
1062 * Y - DOUBLE PRECISION array of dimension at least
1063 * ( 1 + ( n - 1 )*abs( INCY ) ).
1064 * Before entry, the incremented array Y must contain the n
1065 * element vector y. On exit, Y is overwritten by the updated
1066 * vector y.
1067 *
1068 * INCY - INTEGER.
1069 * On entry, INCY specifies the increment for the elements of
1070 * Y. INCY must not be zero.
1071 * Unchanged on exit.
1072 *
1073 *
1074 * Level 2 Blas routine.
1075 *
1076 * -- Written on 22-October-1986.
1077 * Jack Dongarra, Argonne National Lab.
1078 * Jeremy Du Croz, Nag Central Office.
1079 * Sven Hammarling, Nag Central Office.
1080 * Richard Hanson, Sandia National Labs.
1081 *
1082 *
1083 * .. Parameters ..
1084 **/
C_DSPMV(char uplo,int n,double alpha,double * ap,double * x,int incx,double beta,double * y,int incy)1085 void C_DSPMV(char uplo, int n, double alpha, double *ap, double *x, int incx,
1086 double beta, double *y, int incy)
1087 {
1088 if (n == 0)
1089 return;
1090 if (uplo == 'U' || uplo == 'u')
1091 uplo = 'L';
1092 else if (uplo == 'L' || uplo == 'l')
1093 uplo = 'U';
1094 else
1095 throw std::invalid_argument("C_DSPMV uplo argument is invalid.");
1096 ::F_DSPMV(&uplo, &n, &alpha, ap, x, &incx, &beta, y, &incy);
1097 }
1098
1099 /**
1100 * Purpose
1101 * =======
1102 *
1103 * DSPR performs the symmetric rank 1 operation
1104 *
1105 * A := alpha*x*x' + A,
1106 *
1107 * where alpha is a real scalar, x is an n element vector and A is an
1108 * n by n symmetric matrix, supplied in packed form.
1109 *
1110 * Arguments
1111 * ==========
1112 *
1113 * UPLO - CHARACTER*1.
1114 * On entry, UPLO specifies whether the upper or lower
1115 * triangular part of the matrix A is supplied in the packed
1116 * array AP as follows:
1117 *
1118 * UPLO = 'U' or 'u' The upper triangular part of A is
1119 * supplied in AP.
1120 *
1121 * UPLO = 'L' or 'l' The lower triangular part of A is
1122 * supplied in AP.
1123 *
1124 * Unchanged on exit.
1125 *
1126 * N - INTEGER.
1127 * On entry, N specifies the order of the matrix A.
1128 * N must be at least zero.
1129 * Unchanged on exit.
1130 *
1131 * ALPHA - DOUBLE PRECISION.
1132 * On entry, ALPHA specifies the scalar alpha.
1133 * Unchanged on exit.
1134 *
1135 * X - DOUBLE PRECISION array of dimension at least
1136 * ( 1 + ( n - 1 )*abs( INCX ) ).
1137 * Before entry, the incremented array X must contain the n
1138 * element vector x.
1139 * Unchanged on exit.
1140 *
1141 * INCX - INTEGER.
1142 * On entry, INCX specifies the increment for the elements of
1143 * X. INCX must not be zero.
1144 * Unchanged on exit.
1145 *
1146 * AP - DOUBLE PRECISION array of DIMENSION at least
1147 * ( ( n*( n + 1 ) )/2 ).
1148 * Before entry with UPLO = 'U' or 'u', the array AP must
1149 * contain the upper triangular part of the symmetric matrix
1150 * packed sequentially, column by column, so that AP( 1 )
1151 * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
1152 * and a( 2, 2 ) respectively, and so on. On exit, the array
1153 * AP is overwritten by the upper triangular part of the
1154 * updated matrix.
1155 * Before entry with UPLO = 'L' or 'l', the array AP must
1156 * contain the lower triangular part of the symmetric matrix
1157 * packed sequentially, column by column, so that AP( 1 )
1158 * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
1159 * and a( 3, 1 ) respectively, and so on. On exit, the array
1160 * AP is overwritten by the lower triangular part of the
1161 * updated matrix.
1162 *
1163 *
1164 * Level 2 Blas routine.
1165 *
1166 * -- Written on 22-October-1986.
1167 * Jack Dongarra, Argonne National Lab.
1168 * Jeremy Du Croz, Nag Central Office.
1169 * Sven Hammarling, Nag Central Office.
1170 * Richard Hanson, Sandia National Labs.
1171 *
1172 *
1173 * .. Parameters ..
1174 **/
C_DSPR(char uplo,int n,double alpha,double * x,int incx,double * ap)1175 void C_DSPR(char uplo, int n, double alpha, double *x, int incx, double *ap)
1176 {
1177 if (n == 0)
1178 return;
1179 if (uplo == 'U' || uplo == 'u')
1180 uplo = 'L';
1181 else if (uplo == 'L' || uplo == 'l')
1182 uplo = 'U';
1183 else
1184 throw std::invalid_argument("C_DSPR uplo argument is invalid.");
1185 ::F_DSPR(&uplo, &n, &alpha, x, &incx, ap);
1186 }
1187
1188 /**
1189 * Purpose
1190 * =======
1191 *
1192 * DSPR2 performs the symmetric rank 2 operation
1193 *
1194 * A := alpha*x*y' + alpha*y*x' + A,
1195 *
1196 * where alpha is a scalar, x and y are n element vectors and A is an
1197 * n by n symmetric matrix, supplied in packed form.
1198 *
1199 * Arguments
1200 * ==========
1201 *
1202 * UPLO - CHARACTER*1.
1203 * On entry, UPLO specifies whether the upper or lower
1204 * triangular part of the matrix A is supplied in the packed
1205 * array AP as follows:
1206 *
1207 * UPLO = 'U' or 'u' The upper triangular part of A is
1208 * supplied in AP.
1209 *
1210 * UPLO = 'L' or 'l' The lower triangular part of A is
1211 * supplied in AP.
1212 *
1213 * Unchanged on exit.
1214 *
1215 * N - INTEGER.
1216 * On entry, N specifies the order of the matrix A.
1217 * N must be at least zero.
1218 * Unchanged on exit.
1219 *
1220 * ALPHA - DOUBLE PRECISION.
1221 * On entry, ALPHA specifies the scalar alpha.
1222 * Unchanged on exit.
1223 *
1224 * X - DOUBLE PRECISION array of dimension at least
1225 * ( 1 + ( n - 1 )*abs( INCX ) ).
1226 * Before entry, the incremented array X must contain the n
1227 * element vector x.
1228 * Unchanged on exit.
1229 *
1230 * INCX - INTEGER.
1231 * On entry, INCX specifies the increment for the elements of
1232 * X. INCX must not be zero.
1233 * Unchanged on exit.
1234 *
1235 * Y - DOUBLE PRECISION array of dimension at least
1236 * ( 1 + ( n - 1 )*abs( INCY ) ).
1237 * Before entry, the incremented array Y must contain the n
1238 * element vector y.
1239 * Unchanged on exit.
1240 *
1241 * INCY - INTEGER.
1242 * On entry, INCY specifies the increment for the elements of
1243 * Y. INCY must not be zero.
1244 * Unchanged on exit.
1245 *
1246 * AP - DOUBLE PRECISION array of DIMENSION at least
1247 * ( ( n*( n + 1 ) )/2 ).
1248 * Before entry with UPLO = 'U' or 'u', the array AP must
1249 * contain the upper triangular part of the symmetric matrix
1250 * packed sequentially, column by column, so that AP( 1 )
1251 * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
1252 * and a( 2, 2 ) respectively, and so on. On exit, the array
1253 * AP is overwritten by the upper triangular part of the
1254 * updated matrix.
1255 * Before entry with UPLO = 'L' or 'l', the array AP must
1256 * contain the lower triangular part of the symmetric matrix
1257 * packed sequentially, column by column, so that AP( 1 )
1258 * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
1259 * and a( 3, 1 ) respectively, and so on. On exit, the array
1260 * AP is overwritten by the lower triangular part of the
1261 * updated matrix.
1262 *
1263 *
1264 * Level 2 Blas routine.
1265 *
1266 * -- Written on 22-October-1986.
1267 * Jack Dongarra, Argonne National Lab.
1268 * Jeremy Du Croz, Nag Central Office.
1269 * Sven Hammarling, Nag Central Office.
1270 * Richard Hanson, Sandia National Labs.
1271 *
1272 *
1273 * .. Parameters ..
1274 **/
C_DSPR2(char uplo,int n,double alpha,double * x,int incx,double * y,int incy,double * ap)1275 void C_DSPR2(char uplo, int n, double alpha, double *x, int incx, double *y,
1276 int incy, double *ap)
1277 {
1278 if (n == 0)
1279 return;
1280 if (uplo == 'U' || uplo == 'u')
1281 uplo = 'L';
1282 else if (uplo == 'L' || uplo == 'l')
1283 uplo = 'U';
1284 else
1285 throw std::invalid_argument("C_DSPR2 uplo argument is invalid.");
1286 ::F_DSPR2(&uplo, &n, &alpha, x, &incx, y, &incy, ap);
1287 }
1288
1289 /**
1290 * Purpose
1291 * =======
1292 *
1293 * DSYMM performs one of the matrix-matrix operations
1294 *
1295 * C := alpha*A*B + beta*C,
1296 *
1297 * or
1298 *
1299 * C := alpha*B*A + beta*C,
1300 *
1301 * where alpha and beta are scalars, A is a symmetric matrix and B and
1302 * C are m by n matrices.
1303 *
1304 * Arguments
1305 * ==========
1306 *
1307 * SIDE - CHARACTER*1.
1308 * On entry, SIDE specifies whether the symmetric matrix A
1309 * appears on the left or right in the operation as follows:
1310 *
1311 * SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
1312 *
1313 * SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
1314 *
1315 * Unchanged on exit.
1316 *
1317 * UPLO - CHARACTER*1.
1318 * On entry, UPLO specifies whether the upper or lower
1319 * triangular part of the symmetric matrix A is to be
1320 * referenced as follows:
1321 *
1322 * UPLO = 'U' or 'u' Only the upper triangular part of the
1323 * symmetric matrix is to be referenced.
1324 *
1325 * UPLO = 'L' or 'l' Only the lower triangular part of the
1326 * symmetric matrix is to be referenced.
1327 *
1328 * Unchanged on exit.
1329 *
1330 * M - INTEGER.
1331 * On entry, M specifies the number of rows of the matrix C.
1332 * M must be at least zero.
1333 * Unchanged on exit.
1334 *
1335 * N - INTEGER.
1336 * On entry, N specifies the number of columns of the matrix C.
1337 * N must be at least zero.
1338 * Unchanged on exit.
1339 *
1340 * ALPHA - DOUBLE PRECISION.
1341 * On entry, ALPHA specifies the scalar alpha.
1342 * Unchanged on exit.
1343 *
1344 * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
1345 * m when SIDE = 'L' or 'l' and is n otherwise.
1346 * Before entry with SIDE = 'L' or 'l', the m by m part of
1347 * the array A must contain the symmetric matrix, such that
1348 * when UPLO = 'U' or 'u', the leading m by m upper triangular
1349 * part of the array A must contain the upper triangular part
1350 * of the symmetric matrix and the strictly lower triangular
1351 * part of A is not referenced, and when UPLO = 'L' or 'l',
1352 * the leading m by m lower triangular part of the array A
1353 * must contain the lower triangular part of the symmetric
1354 * matrix and the strictly upper triangular part of A is not
1355 * referenced.
1356 * Before entry with SIDE = 'R' or 'r', the n by n part of
1357 * the array A must contain the symmetric matrix, such that
1358 * when UPLO = 'U' or 'u', the leading n by n upper triangular
1359 * part of the array A must contain the upper triangular part
1360 * of the symmetric matrix and the strictly lower triangular
1361 * part of A is not referenced, and when UPLO = 'L' or 'l',
1362 * the leading n by n lower triangular part of the array A
1363 * must contain the lower triangular part of the symmetric
1364 * matrix and the strictly upper triangular part of A is not
1365 * referenced.
1366 * Unchanged on exit.
1367 *
1368 * LDA - INTEGER.
1369 * On entry, LDA specifies the first dimension of A as declared
1370 * in the calling (sub) program. When SIDE = 'L' or 'l' then
1371 * LDA must be at least max( 1, m ), otherwise LDA must be at
1372 * least max( 1, n ).
1373 * Unchanged on exit.
1374 *
1375 * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
1376 * Before entry, the leading m by n part of the array B must
1377 * contain the matrix B.
1378 * Unchanged on exit.
1379 *
1380 * LDB - INTEGER.
1381 * On entry, LDB specifies the first dimension of B as declared
1382 * in the calling (sub) program. LDB must be at least
1383 * max( 1, m ).
1384 * Unchanged on exit.
1385 *
1386 * BETA - DOUBLE PRECISION.
1387 * On entry, BETA specifies the scalar beta. When BETA is
1388 * supplied as zero then C need not be set on input.
1389 * Unchanged on exit.
1390 *
1391 * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
1392 * Before entry, the leading m by n part of the array C must
1393 * contain the matrix C, except when beta is zero, in which
1394 * case C need not be set on entry.
1395 * On exit, the array C is overwritten by the m by n updated
1396 * matrix.
1397 *
1398 * LDC - INTEGER.
1399 * On entry, LDC specifies the first dimension of C as declared
1400 * in the calling (sub) program. LDC must be at least
1401 * max( 1, m ).
1402 * Unchanged on exit.
1403 *
1404 *
1405 * Level 3 Blas routine.
1406 *
1407 * -- Written on 8-February-1989.
1408 * Jack Dongarra, Argonne National Laboratory.
1409 * Iain Duff, AERE Harwell.
1410 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
1411 * Sven Hammarling, Numerical Algorithms Group Ltd.
1412 *
1413 *
1414 * .. External Functions ..
1415 **/
C_DSYMM(char side,char uplo,int m,int n,double alpha,double * a,int lda,double * b,int ldb,double beta,double * c,int ldc)1416 void C_DSYMM(char side, char uplo, int m, int n, double alpha, double *a,
1417 int lda, double *b, int ldb, double beta, double *c, int ldc)
1418 {
1419 if (m == 0 || n == 0)
1420 return;
1421 if (uplo == 'U' || uplo == 'u')
1422 uplo = 'L';
1423 else if (uplo == 'L' || uplo == 'l')
1424 uplo = 'U';
1425 else
1426 throw std::invalid_argument("C_DSYMM uplo argument is invalid.");
1427 if (side == 'L' || side == 'L')
1428 side = 'R';
1429 else if (side == 'R' || side == 'r')
1430 side = 'L';
1431 else
1432 throw std::invalid_argument("C_DSYMM side argument is invalid.");
1433 ::F_DSYMM(&side, &uplo, &n, &m, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
1434 }
1435
1436 /**
1437 * Purpose
1438 * =======
1439 *
1440 * DSYMV performs the matrix-vector operation
1441 *
1442 * y := alpha*A*x + beta*y,
1443 *
1444 * where alpha and beta are scalars, x and y are n element vectors and
1445 * A is an n by n symmetric matrix.
1446 *
1447 * Arguments
1448 * ==========
1449 *
1450 * UPLO - CHARACTER*1.
1451 * On entry, UPLO specifies whether the upper or lower
1452 * triangular part of the array A is to be referenced as
1453 * follows:
1454 *
1455 * UPLO = 'U' or 'u' Only the upper triangular part of A
1456 * is to be referenced.
1457 *
1458 * UPLO = 'L' or 'l' Only the lower triangular part of A
1459 * is to be referenced.
1460 *
1461 * Unchanged on exit.
1462 *
1463 * N - INTEGER.
1464 * On entry, N specifies the order of the matrix A.
1465 * N must be at least zero.
1466 * Unchanged on exit.
1467 *
1468 * ALPHA - DOUBLE PRECISION.
1469 * On entry, ALPHA specifies the scalar alpha.
1470 * Unchanged on exit.
1471 *
1472 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
1473 * Before entry with UPLO = 'U' or 'u', the leading n by n
1474 * upper triangular part of the array A must contain the upper
1475 * triangular part of the symmetric matrix and the strictly
1476 * lower triangular part of A is not referenced.
1477 * Before entry with UPLO = 'L' or 'l', the leading n by n
1478 * lower triangular part of the array A must contain the lower
1479 * triangular part of the symmetric matrix and the strictly
1480 * upper triangular part of A is not referenced.
1481 * Unchanged on exit.
1482 *
1483 * LDA - INTEGER.
1484 * On entry, LDA specifies the first dimension of A as declared
1485 * in the calling (sub) program. LDA must be at least
1486 * max( 1, n ).
1487 * Unchanged on exit.
1488 *
1489 * X - DOUBLE PRECISION array of dimension at least
1490 * ( 1 + ( n - 1 )*abs( INCX ) ).
1491 * Before entry, the incremented array X must contain the n
1492 * element vector x.
1493 * Unchanged on exit.
1494 *
1495 * INCX - INTEGER.
1496 * On entry, INCX specifies the increment for the elements of
1497 * X. INCX must not be zero.
1498 * Unchanged on exit.
1499 *
1500 * BETA - DOUBLE PRECISION.
1501 * On entry, BETA specifies the scalar beta. When BETA is
1502 * supplied as zero then Y need not be set on input.
1503 * Unchanged on exit.
1504 *
1505 * Y - DOUBLE PRECISION array of dimension at least
1506 * ( 1 + ( n - 1 )*abs( INCY ) ).
1507 * Before entry, the incremented array Y must contain the n
1508 * element vector y. On exit, Y is overwritten by the updated
1509 * vector y.
1510 *
1511 * INCY - INTEGER.
1512 * On entry, INCY specifies the increment for the elements of
1513 * Y. INCY must not be zero.
1514 * Unchanged on exit.
1515 *
1516 *
1517 * Level 2 Blas routine.
1518 *
1519 * -- Written on 22-October-1986.
1520 * Jack Dongarra, Argonne National Lab.
1521 * Jeremy Du Croz, Nag Central Office.
1522 * Sven Hammarling, Nag Central Office.
1523 * Richard Hanson, Sandia National Labs.
1524 *
1525 *
1526 * .. Parameters ..
1527 **/
C_DSYMV(char uplo,int n,double alpha,double * a,int lda,double * x,int incx,double beta,double * y,int incy)1528 void C_DSYMV(char uplo, int n, double alpha, double *a, int lda, double *x,
1529 int incx, double beta, double *y, int incy)
1530 {
1531 if (n == 0)
1532 return;
1533 if (uplo == 'U' || uplo == 'u')
1534 uplo = 'L';
1535 else if (uplo == 'L' || uplo == 'l')
1536 uplo = 'U';
1537 else
1538 throw std::invalid_argument("C_DSYMV uplo argument is invalid.");
1539 ::F_DSYMV(&uplo, &n, &alpha, a, &lda, x, &incx, &beta, y, &incy);
1540 }
1541
1542 /**
1543 * Purpose
1544 * =======
1545 *
1546 * DSYR performs the symmetric rank 1 operation
1547 *
1548 * A := alpha*x*x' + A,
1549 *
1550 * where alpha is a real scalar, x is an n element vector and A is an
1551 * n by n symmetric matrix.
1552 *
1553 * Arguments
1554 * ==========
1555 *
1556 * UPLO - CHARACTER*1.
1557 * On entry, UPLO specifies whether the upper or lower
1558 * triangular part of the array A is to be referenced as
1559 * follows:
1560 *
1561 * UPLO = 'U' or 'u' Only the upper triangular part of A
1562 * is to be referenced.
1563 *
1564 * UPLO = 'L' or 'l' Only the lower triangular part of A
1565 * is to be referenced.
1566 *
1567 * Unchanged on exit.
1568 *
1569 * N - INTEGER.
1570 * On entry, N specifies the order of the matrix A.
1571 * N must be at least zero.
1572 * Unchanged on exit.
1573 *
1574 * ALPHA - DOUBLE PRECISION.
1575 * On entry, ALPHA specifies the scalar alpha.
1576 * Unchanged on exit.
1577 *
1578 * X - DOUBLE PRECISION array of dimension at least
1579 * ( 1 + ( n - 1 )*abs( INCX ) ).
1580 * Before entry, the incremented array X must contain the n
1581 * element vector x.
1582 * Unchanged on exit.
1583 *
1584 * INCX - INTEGER.
1585 * On entry, INCX specifies the increment for the elements of
1586 * X. INCX must not be zero.
1587 * Unchanged on exit.
1588 *
1589 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
1590 * Before entry with UPLO = 'U' or 'u', the leading n by n
1591 * upper triangular part of the array A must contain the upper
1592 * triangular part of the symmetric matrix and the strictly
1593 * lower triangular part of A is not referenced. On exit, the
1594 * upper triangular part of the array A is overwritten by the
1595 * upper triangular part of the updated matrix.
1596 * Before entry with UPLO = 'L' or 'l', the leading n by n
1597 * lower triangular part of the array A must contain the lower
1598 * triangular part of the symmetric matrix and the strictly
1599 * upper triangular part of A is not referenced. On exit, the
1600 * lower triangular part of the array A is overwritten by the
1601 * lower triangular part of the updated matrix.
1602 *
1603 * LDA - INTEGER.
1604 * On entry, LDA specifies the first dimension of A as declared
1605 * in the calling (sub) program. LDA must be at least
1606 * max( 1, n ).
1607 * Unchanged on exit.
1608 *
1609 *
1610 * Level 2 Blas routine.
1611 *
1612 * -- Written on 22-October-1986.
1613 * Jack Dongarra, Argonne National Lab.
1614 * Jeremy Du Croz, Nag Central Office.
1615 * Sven Hammarling, Nag Central Office.
1616 * Richard Hanson, Sandia National Labs.
1617 *
1618 *
1619 * .. Parameters ..
1620 **/
C_DSYR(char uplo,int n,double alpha,double * x,int incx,double * a,int lda)1621 void C_DSYR(char uplo, int n, double alpha, double *x, int incx, double *a,
1622 int lda)
1623 {
1624 if (n == 0)
1625 return;
1626 if (uplo == 'U' || uplo == 'u')
1627 uplo = 'L';
1628 else if (uplo == 'L' || uplo == 'l')
1629 uplo = 'U';
1630 else
1631 throw std::invalid_argument("C_DSYR uplo argument is invalid.");
1632 ::F_DSYR(&uplo, &n, &alpha, x, &incx, a, &lda);
1633 }
1634
1635 /**
1636 * Purpose
1637 * =======
1638 *
1639 * DSYR2 performs the symmetric rank 2 operation
1640 *
1641 * A := alpha*x*y' + alpha*y*x' + A,
1642 *
1643 * where alpha is a scalar, x and y are n element vectors and A is an n
1644 * by n symmetric matrix.
1645 *
1646 * Arguments
1647 * ==========
1648 *
1649 * UPLO - CHARACTER*1.
1650 * On entry, UPLO specifies whether the upper or lower
1651 * triangular part of the array A is to be referenced as
1652 * follows:
1653 *
1654 * UPLO = 'U' or 'u' Only the upper triangular part of A
1655 * is to be referenced.
1656 *
1657 * UPLO = 'L' or 'l' Only the lower triangular part of A
1658 * is to be referenced.
1659 *
1660 * Unchanged on exit.
1661 *
1662 * N - INTEGER.
1663 * On entry, N specifies the order of the matrix A.
1664 * N must be at least zero.
1665 * Unchanged on exit.
1666 *
1667 * ALPHA - DOUBLE PRECISION.
1668 * On entry, ALPHA specifies the scalar alpha.
1669 * Unchanged on exit.
1670 *
1671 * X - DOUBLE PRECISION array of dimension at least
1672 * ( 1 + ( n - 1 )*abs( INCX ) ).
1673 * Before entry, the incremented array X must contain the n
1674 * element vector x.
1675 * Unchanged on exit.
1676 *
1677 * INCX - INTEGER.
1678 * On entry, INCX specifies the increment for the elements of
1679 * X. INCX must not be zero.
1680 * Unchanged on exit.
1681 *
1682 * Y - DOUBLE PRECISION array of dimension at least
1683 * ( 1 + ( n - 1 )*abs( INCY ) ).
1684 * Before entry, the incremented array Y must contain the n
1685 * element vector y.
1686 * Unchanged on exit.
1687 *
1688 * INCY - INTEGER.
1689 * On entry, INCY specifies the increment for the elements of
1690 * Y. INCY must not be zero.
1691 * Unchanged on exit.
1692 *
1693 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
1694 * Before entry with UPLO = 'U' or 'u', the leading n by n
1695 * upper triangular part of the array A must contain the upper
1696 * triangular part of the symmetric matrix and the strictly
1697 * lower triangular part of A is not referenced. On exit, the
1698 * upper triangular part of the array A is overwritten by the
1699 * upper triangular part of the updated matrix.
1700 * Before entry with UPLO = 'L' or 'l', the leading n by n
1701 * lower triangular part of the array A must contain the lower
1702 * triangular part of the symmetric matrix and the strictly
1703 * upper triangular part of A is not referenced. On exit, the
1704 * lower triangular part of the array A is overwritten by the
1705 * lower triangular part of the updated matrix.
1706 *
1707 * LDA - INTEGER.
1708 * On entry, LDA specifies the first dimension of A as declared
1709 * in the calling (sub) program. LDA must be at least
1710 * max( 1, n ).
1711 * Unchanged on exit.
1712 *
1713 *
1714 * Level 2 Blas routine.
1715 *
1716 * -- Written on 22-October-1986.
1717 * Jack Dongarra, Argonne National Lab.
1718 * Jeremy Du Croz, Nag Central Office.
1719 * Sven Hammarling, Nag Central Office.
1720 * Richard Hanson, Sandia National Labs.
1721 *
1722 *
1723 * .. Parameters ..
1724 **/
C_DSYR2(char uplo,int n,double alpha,double * x,int incx,double * y,int incy,double * a,int lda)1725 void C_DSYR2(char uplo, int n, double alpha, double *x, int incx, double *y,
1726 int incy, double *a, int lda)
1727 {
1728 if (n == 0)
1729 return;
1730 if (uplo == 'U' || uplo == 'u')
1731 uplo = 'L';
1732 else if (uplo == 'L' || uplo == 'l')
1733 uplo = 'U';
1734 else
1735 throw std::invalid_argument("C_DSYR2 uplo argument is invalid.");
1736 ::F_DSYR2(&uplo, &n, &alpha, x, &incx, y, &incy, a, &lda);
1737 }
1738
1739 /**
1740 * Purpose
1741 * =======
1742 *
1743 * DSYR2K performs one of the symmetric rank 2k operations
1744 *
1745 * C := alpha*A*B' + alpha*B*A' + beta*C,
1746 *
1747 * or
1748 *
1749 * C := alpha*A'*B + alpha*B'*A + beta*C,
1750 *
1751 * where alpha and beta are scalars, C is an n by n symmetric matrix
1752 * and A and B are n by k matrices in the first case and k by n
1753 * matrices in the second case.
1754 *
1755 * Arguments
1756 * ==========
1757 *
1758 * UPLO - CHARACTER*1.
1759 * On entry, UPLO specifies whether the upper or lower
1760 * triangular part of the array C is to be referenced as
1761 * follows:
1762 *
1763 * UPLO = 'U' or 'u' Only the upper triangular part of C
1764 * is to be referenced.
1765 *
1766 * UPLO = 'L' or 'l' Only the lower triangular part of C
1767 * is to be referenced.
1768 *
1769 * Unchanged on exit.
1770 *
1771 * TRANS - CHARACTER*1.
1772 * On entry, TRANS specifies the operation to be performed as
1773 * follows:
1774 *
1775 * TRANS = 'N' or 'n' C := alpha*A*B' + alpha*B*A' +
1776 * beta*C.
1777 *
1778 * TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*A +
1779 * beta*C.
1780 *
1781 * TRANS = 'C' or 'c' C := alpha*A'*B + alpha*B'*A +
1782 * beta*C.
1783 *
1784 * Unchanged on exit.
1785 *
1786 * N - INTEGER.
1787 * On entry, N specifies the order of the matrix C. N must be
1788 * at least zero.
1789 * Unchanged on exit.
1790 *
1791 * K - INTEGER.
1792 * On entry with TRANS = 'N' or 'n', K specifies the number
1793 * of columns of the matrices A and B, and on entry with
1794 * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
1795 * of rows of the matrices A and B. K must be at least zero.
1796 * Unchanged on exit.
1797 *
1798 * ALPHA - DOUBLE PRECISION.
1799 * On entry, ALPHA specifies the scalar alpha.
1800 * Unchanged on exit.
1801 *
1802 * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
1803 * k when TRANS = 'N' or 'n', and is n otherwise.
1804 * Before entry with TRANS = 'N' or 'n', the leading n by k
1805 * part of the array A must contain the matrix A, otherwise
1806 * the leading k by n part of the array A must contain the
1807 * matrix A.
1808 * Unchanged on exit.
1809 *
1810 * LDA - INTEGER.
1811 * On entry, LDA specifies the first dimension of A as declared
1812 * in the calling (sub) program. When TRANS = 'N' or 'n'
1813 * then LDA must be at least max( 1, n ), otherwise LDA must
1814 * be at least max( 1, k ).
1815 * Unchanged on exit.
1816 *
1817 * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
1818 * k when TRANS = 'N' or 'n', and is n otherwise.
1819 * Before entry with TRANS = 'N' or 'n', the leading n by k
1820 * part of the array B must contain the matrix B, otherwise
1821 * the leading k by n part of the array B must contain the
1822 * matrix B.
1823 * Unchanged on exit.
1824 *
1825 * LDB - INTEGER.
1826 * On entry, LDB specifies the first dimension of B as declared
1827 * in the calling (sub) program. When TRANS = 'N' or 'n'
1828 * then LDB must be at least max( 1, n ), otherwise LDB must
1829 * be at least max( 1, k ).
1830 * Unchanged on exit.
1831 *
1832 * BETA - DOUBLE PRECISION.
1833 * On entry, BETA specifies the scalar beta.
1834 * Unchanged on exit.
1835 *
1836 * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
1837 * Before entry with UPLO = 'U' or 'u', the leading n by n
1838 * upper triangular part of the array C must contain the upper
1839 * triangular part of the symmetric matrix and the strictly
1840 * lower triangular part of C is not referenced. On exit, the
1841 * upper triangular part of the array C is overwritten by the
1842 * upper triangular part of the updated matrix.
1843 * Before entry with UPLO = 'L' or 'l', the leading n by n
1844 * lower triangular part of the array C must contain the lower
1845 * triangular part of the symmetric matrix and the strictly
1846 * upper triangular part of C is not referenced. On exit, the
1847 * lower triangular part of the array C is overwritten by the
1848 * lower triangular part of the updated matrix.
1849 *
1850 * LDC - INTEGER.
1851 * On entry, LDC specifies the first dimension of C as declared
1852 * in the calling (sub) program. LDC must be at least
1853 * max( 1, n ).
1854 * Unchanged on exit.
1855 *
1856 *
1857 * Level 3 Blas routine.
1858 *
1859 *
1860 * -- Written on 8-February-1989.
1861 * Jack Dongarra, Argonne National Laboratory.
1862 * Iain Duff, AERE Harwell.
1863 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
1864 * Sven Hammarling, Numerical Algorithms Group Ltd.
1865 *
1866 *
1867 * .. External Functions ..
1868 **/
C_DSYR2K(char uplo,char trans,int n,int k,double alpha,double * a,int lda,double * b,int ldb,double beta,double * c,int ldc)1869 void C_DSYR2K(char uplo, char trans, int n, int k, double alpha, double *a,
1870 int lda, double *b, int ldb, double beta, double *c, int ldc)
1871 {
1872 if (n == 0 || k == 0)
1873 return;
1874 if (uplo == 'U' || uplo == 'u')
1875 uplo = 'L';
1876 else if (uplo == 'L' || uplo == 'l')
1877 uplo = 'U';
1878 else
1879 throw std::invalid_argument("C_DSYR2K uplo argument is invalid.");
1880 if (trans == 'N' || trans == 'n')
1881 trans = 'T';
1882 else if (trans == 'T' || trans == 't')
1883 trans = 'N';
1884 else
1885 throw std::invalid_argument("C_DSYR2K trans argument is invalid.");
1886
1887 ::F_DSYR2K(&uplo, &trans, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
1888 }
1889
1890 /**
1891 * Purpose
1892 * =======
1893 *
1894 * DSYRK performs one of the symmetric rank k operations
1895 *
1896 * C := alpha*A*A' + beta*C,
1897 *
1898 * or
1899 *
1900 * C := alpha*A'*A + beta*C,
1901 *
1902 * where alpha and beta are scalars, C is an n by n symmetric matrix
1903 * and A is an n by k matrix in the first case and a k by n matrix
1904 * in the second case.
1905 *
1906 * Arguments
1907 * ==========
1908 *
1909 * UPLO - CHARACTER*1.
1910 * On entry, UPLO specifies whether the upper or lower
1911 * triangular part of the array C is to be referenced as
1912 * follows:
1913 *
1914 * UPLO = 'U' or 'u' Only the upper triangular part of C
1915 * is to be referenced.
1916 *
1917 * UPLO = 'L' or 'l' Only the lower triangular part of C
1918 * is to be referenced.
1919 *
1920 * Unchanged on exit.
1921 *
1922 * TRANS - CHARACTER*1.
1923 * On entry, TRANS specifies the operation to be performed as
1924 * follows:
1925 *
1926 * TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
1927 *
1928 * TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
1929 *
1930 * TRANS = 'C' or 'c' C := alpha*A'*A + beta*C.
1931 *
1932 * Unchanged on exit.
1933 *
1934 * N - INTEGER.
1935 * On entry, N specifies the order of the matrix C. N must be
1936 * at least zero.
1937 * Unchanged on exit.
1938 *
1939 * K - INTEGER.
1940 * On entry with TRANS = 'N' or 'n', K specifies the number
1941 * of columns of the matrix A, and on entry with
1942 * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
1943 * of rows of the matrix A. K must be at least zero.
1944 * Unchanged on exit.
1945 *
1946 * ALPHA - DOUBLE PRECISION.
1947 * On entry, ALPHA specifies the scalar alpha.
1948 * Unchanged on exit.
1949 *
1950 * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
1951 * k when TRANS = 'N' or 'n', and is n otherwise.
1952 * Before entry with TRANS = 'N' or 'n', the leading n by k
1953 * part of the array A must contain the matrix A, otherwise
1954 * the leading k by n part of the array A must contain the
1955 * matrix A.
1956 * Unchanged on exit.
1957 *
1958 * LDA - INTEGER.
1959 * On entry, LDA specifies the first dimension of A as declared
1960 * in the calling (sub) program. When TRANS = 'N' or 'n'
1961 * then LDA must be at least max( 1, n ), otherwise LDA must
1962 * be at least max( 1, k ).
1963 * Unchanged on exit.
1964 *
1965 * BETA - DOUBLE PRECISION.
1966 * On entry, BETA specifies the scalar beta.
1967 * Unchanged on exit.
1968 *
1969 * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
1970 * Before entry with UPLO = 'U' or 'u', the leading n by n
1971 * upper triangular part of the array C must contain the upper
1972 * triangular part of the symmetric matrix and the strictly
1973 * lower triangular part of C is not referenced. On exit, the
1974 * upper triangular part of the array C is overwritten by the
1975 * upper triangular part of the updated matrix.
1976 * Before entry with UPLO = 'L' or 'l', the leading n by n
1977 * lower triangular part of the array C must contain the lower
1978 * triangular part of the symmetric matrix and the strictly
1979 * upper triangular part of C is not referenced. On exit, the
1980 * lower triangular part of the array C is overwritten by the
1981 * lower triangular part of the updated matrix.
1982 *
1983 * LDC - INTEGER.
1984 * On entry, LDC specifies the first dimension of C as declared
1985 * in the calling (sub) program. LDC must be at least
1986 * max( 1, n ).
1987 * Unchanged on exit.
1988 *
1989 *
1990 * Level 3 Blas routine.
1991 *
1992 * -- Written on 8-February-1989.
1993 * Jack Dongarra, Argonne National Laboratory.
1994 * Iain Duff, AERE Harwell.
1995 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
1996 * Sven Hammarling, Numerical Algorithms Group Ltd.
1997 *
1998 *
1999 * .. External Functions ..
2000 **/
C_DSYRK(char uplo,char trans,int n,int k,double alpha,double * a,int lda,double beta,double * c,int ldc)2001 void C_DSYRK(char uplo, char trans, int n, int k, double alpha, double *a,
2002 int lda, double beta, double *c, int ldc)
2003 {
2004 if (n == 0 || k == 0)
2005 return;
2006 if (uplo == 'U' || uplo == 'u')
2007 uplo = 'L';
2008 else if (uplo == 'L' || uplo == 'l')
2009 uplo = 'U';
2010 else
2011 throw std::invalid_argument("C_DSYRK uplo argument is invalid.");
2012 if (trans == 'N' || trans == 'n')
2013 trans = 'T';
2014 else if (trans == 'T' || trans == 't')
2015 trans = 'N';
2016 else
2017 throw std::invalid_argument("C_DSYRK trans argument is invalid.");
2018
2019 ::F_DSYRK(&uplo, &trans, &n, &k, &alpha, a, &lda, &beta, c, &ldc);
2020 }
2021
2022 /**
2023 * Purpose
2024 * =======
2025 *
2026 * DTBMV performs one of the matrix-vector operations
2027 *
2028 * x := A*x, or x := A'*x,
2029 *
2030 * where x is an n element vector and A is an n by n unit, or non-unit,
2031 * upper or lower triangular band matrix, with ( k + 1 ) diagonals.
2032 *
2033 * Arguments
2034 * ==========
2035 *
2036 * UPLO - CHARACTER*1.
2037 * On entry, UPLO specifies whether the matrix is an upper or
2038 * lower triangular matrix as follows:
2039 *
2040 * UPLO = 'U' or 'u' A is an upper triangular matrix.
2041 *
2042 * UPLO = 'L' or 'l' A is a lower triangular matrix.
2043 *
2044 * Unchanged on exit.
2045 *
2046 * TRANS - CHARACTER*1.
2047 * On entry, TRANS specifies the operation to be performed as
2048 * follows:
2049 *
2050 * TRANS = 'N' or 'n' x := A*x.
2051 *
2052 * TRANS = 'T' or 't' x := A'*x.
2053 *
2054 * TRANS = 'C' or 'c' x := A'*x.
2055 *
2056 * Unchanged on exit.
2057 *
2058 * DIAG - CHARACTER*1.
2059 * On entry, DIAG specifies whether or not A is unit
2060 * triangular as follows:
2061 *
2062 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
2063 *
2064 * DIAG = 'N' or 'n' A is not assumed to be unit
2065 * triangular.
2066 *
2067 * Unchanged on exit.
2068 *
2069 * N - INTEGER.
2070 * On entry, N specifies the order of the matrix A.
2071 * N must be at least zero.
2072 * Unchanged on exit.
2073 *
2074 * K - INTEGER.
2075 * On entry with UPLO = 'U' or 'u', K specifies the number of
2076 * super-diagonals of the matrix A.
2077 * On entry with UPLO = 'L' or 'l', K specifies the number of
2078 * sub-diagonals of the matrix A.
2079 * K must satisfy 0 .le. K.
2080 * Unchanged on exit.
2081 *
2082 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
2083 * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
2084 * by n part of the array A must contain the upper triangular
2085 * band part of the matrix of coefficients, supplied column by
2086 * column, with the leading diagonal of the matrix in row
2087 * ( k + 1 ) of the array, the first super-diagonal starting at
2088 * position 2 in row k, and so on. The top left k by k triangle
2089 * of the array A is not referenced.
2090 * The following program segment will transfer an upper
2091 * triangular band matrix from conventional full matrix storage
2092 * to band storage:
2093 *
2094 * DO 20, J = 1, N
2095 * M = K + 1 - J
2096 * DO 10, I = MAX( 1, J - K ), J
2097 * A( M + I, J ) = matrix( I, J )
2098 * 10 CONTINUE
2099 * 20 CONTINUE
2100 *
2101 * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
2102 * by n part of the array A must contain the lower triangular
2103 * band part of the matrix of coefficients, supplied column by
2104 * column, with the leading diagonal of the matrix in row 1 of
2105 * the array, the first sub-diagonal starting at position 1 in
2106 * row 2, and so on. The bottom right k by k triangle of the
2107 * array A is not referenced.
2108 * The following program segment will transfer a lower
2109 * triangular band matrix from conventional full matrix storage
2110 * to band storage:
2111 *
2112 * DO 20, J = 1, N
2113 * M = 1 - J
2114 * DO 10, I = J, MIN( N, J + K )
2115 * A( M + I, J ) = matrix( I, J )
2116 * 10 CONTINUE
2117 * 20 CONTINUE
2118 *
2119 * Note that when DIAG = 'U' or 'u' the elements of the array A
2120 * corresponding to the diagonal elements of the matrix are not
2121 * referenced, but are assumed to be unity.
2122 * Unchanged on exit.
2123 *
2124 * LDA - INTEGER.
2125 * On entry, LDA specifies the first dimension of A as declared
2126 * in the calling (sub) program. LDA must be at least
2127 * ( k + 1 ).
2128 * Unchanged on exit.
2129 *
2130 * X - DOUBLE PRECISION array of dimension at least
2131 * ( 1 + ( n - 1 )*abs( INCX ) ).
2132 * Before entry, the incremented array X must contain the n
2133 * element vector x. On exit, X is overwritten with the
2134 * tranformed vector x.
2135 *
2136 * INCX - INTEGER.
2137 * On entry, INCX specifies the increment for the elements of
2138 * X. INCX must not be zero.
2139 * Unchanged on exit.
2140 *
2141 *
2142 * Level 2 Blas routine.
2143 *
2144 * -- Written on 22-October-1986.
2145 * Jack Dongarra, Argonne National Lab.
2146 * Jeremy Du Croz, Nag Central Office.
2147 * Sven Hammarling, Nag Central Office.
2148 * Richard Hanson, Sandia National Labs.
2149 *
2150 *
2151 * .. Parameters ..
2152 **/
C_DTBMV(char uplo,char trans,char diag,int n,int k,double * a,int lda,double * x,int incx)2153 void C_DTBMV(char uplo, char trans, char diag, int n, int k, double *a, int lda,
2154 double *x, int incx)
2155 {
2156 if (n == 0)
2157 return;
2158 if (uplo == 'U' || uplo == 'u')
2159 uplo = 'L';
2160 else if (uplo == 'L' || uplo == 'l')
2161 uplo = 'U';
2162 else
2163 throw std::invalid_argument("C_DTBMV uplo argument is invalid.");
2164 if (trans == 'N' || trans == 'n')
2165 trans = 'T';
2166 else if (trans == 'T' || trans == 't')
2167 trans = 'N';
2168 else
2169 throw std::invalid_argument("C_DTBMV trans argument is invalid.");
2170 ::F_DTBMV(&uplo, &trans, &diag, &n, &k, a, &lda, x, &incx);
2171 }
2172
2173 /**
2174 * Purpose
2175 * =======
2176 *
2177 * DTBSV solves one of the systems of equations
2178 *
2179 * A*x = b, or A'*x = b,
2180 *
2181 * where b and x are n element vectors and A is an n by n unit, or
2182 * non-unit, upper or lower triangular band matrix, with ( k + 1 )
2183 * diagonals.
2184 *
2185 * No test for singularity or near-singularity is included in this
2186 * routine. Such tests must be performed before calling this routine.
2187 *
2188 * Arguments
2189 * ==========
2190 *
2191 * UPLO - CHARACTER*1.
2192 * On entry, UPLO specifies whether the matrix is an upper or
2193 * lower triangular matrix as follows:
2194 *
2195 * UPLO = 'U' or 'u' A is an upper triangular matrix.
2196 *
2197 * UPLO = 'L' or 'l' A is a lower triangular matrix.
2198 *
2199 * Unchanged on exit.
2200 *
2201 * TRANS - CHARACTER*1.
2202 * On entry, TRANS specifies the equations to be solved as
2203 * follows:
2204 *
2205 * TRANS = 'N' or 'n' A*x = b.
2206 *
2207 * TRANS = 'T' or 't' A'*x = b.
2208 *
2209 * TRANS = 'C' or 'c' A'*x = b.
2210 *
2211 * Unchanged on exit.
2212 *
2213 * DIAG - CHARACTER*1.
2214 * On entry, DIAG specifies whether or not A is unit
2215 * triangular as follows:
2216 *
2217 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
2218 *
2219 * DIAG = 'N' or 'n' A is not assumed to be unit
2220 * triangular.
2221 *
2222 * Unchanged on exit.
2223 *
2224 * N - INTEGER.
2225 * On entry, N specifies the order of the matrix A.
2226 * N must be at least zero.
2227 * Unchanged on exit.
2228 *
2229 * K - INTEGER.
2230 * On entry with UPLO = 'U' or 'u', K specifies the number of
2231 * super-diagonals of the matrix A.
2232 * On entry with UPLO = 'L' or 'l', K specifies the number of
2233 * sub-diagonals of the matrix A.
2234 * K must satisfy 0 .le. K.
2235 * Unchanged on exit.
2236 *
2237 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
2238 * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
2239 * by n part of the array A must contain the upper triangular
2240 * band part of the matrix of coefficients, supplied column by
2241 * column, with the leading diagonal of the matrix in row
2242 * ( k + 1 ) of the array, the first super-diagonal starting at
2243 * position 2 in row k, and so on. The top left k by k triangle
2244 * of the array A is not referenced.
2245 * The following program segment will transfer an upper
2246 * triangular band matrix from conventional full matrix storage
2247 * to band storage:
2248 *
2249 * DO 20, J = 1, N
2250 * M = K + 1 - J
2251 * DO 10, I = MAX( 1, J - K ), J
2252 * A( M + I, J ) = matrix( I, J )
2253 * 10 CONTINUE
2254 * 20 CONTINUE
2255 *
2256 * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
2257 * by n part of the array A must contain the lower triangular
2258 * band part of the matrix of coefficients, supplied column by
2259 * column, with the leading diagonal of the matrix in row 1 of
2260 * the array, the first sub-diagonal starting at position 1 in
2261 * row 2, and so on. The bottom right k by k triangle of the
2262 * array A is not referenced.
2263 * The following program segment will transfer a lower
2264 * triangular band matrix from conventional full matrix storage
2265 * to band storage:
2266 *
2267 * DO 20, J = 1, N
2268 * M = 1 - J
2269 * DO 10, I = J, MIN( N, J + K )
2270 * A( M + I, J ) = matrix( I, J )
2271 * 10 CONTINUE
2272 * 20 CONTINUE
2273 *
2274 * Note that when DIAG = 'U' or 'u' the elements of the array A
2275 * corresponding to the diagonal elements of the matrix are not
2276 * referenced, but are assumed to be unity.
2277 * Unchanged on exit.
2278 *
2279 * LDA - INTEGER.
2280 * On entry, LDA specifies the first dimension of A as declared
2281 * in the calling (sub) program. LDA must be at least
2282 * ( k + 1 ).
2283 * Unchanged on exit.
2284 *
2285 * X - DOUBLE PRECISION array of dimension at least
2286 * ( 1 + ( n - 1 )*abs( INCX ) ).
2287 * Before entry, the incremented array X must contain the n
2288 * element right-hand side vector b. On exit, X is overwritten
2289 * with the solution vector x.
2290 *
2291 * INCX - INTEGER.
2292 * On entry, INCX specifies the increment for the elements of
2293 * X. INCX must not be zero.
2294 * Unchanged on exit.
2295 *
2296 *
2297 * Level 2 Blas routine.
2298 *
2299 * -- Written on 22-October-1986.
2300 * Jack Dongarra, Argonne National Lab.
2301 * Jeremy Du Croz, Nag Central Office.
2302 * Sven Hammarling, Nag Central Office.
2303 * Richard Hanson, Sandia National Labs.
2304 *
2305 *
2306 * .. Parameters ..
2307 **/
C_DTBSV(char uplo,char trans,char diag,int n,int k,double * a,int lda,double * x,int incx)2308 void C_DTBSV(char uplo, char trans, char diag, int n, int k, double *a, int lda,
2309 double *x, int incx)
2310 {
2311 if (n == 0)
2312 return;
2313 if (uplo == 'U' || uplo == 'u')
2314 uplo = 'L';
2315 else if (uplo == 'L' || uplo == 'l')
2316 uplo = 'U';
2317 else
2318 throw std::invalid_argument("C_DTBSV uplo argument is invalid.");
2319 if (trans == 'N' || trans == 'n')
2320 trans = 'T';
2321 else if (trans == 'T' || trans == 't')
2322 trans = 'N';
2323 else
2324 throw std::invalid_argument("C_DTBSV trans argument is invalid.");
2325 ::F_DTBSV(&uplo, &trans, &diag, &n, &k, a, &lda, x, &incx);
2326 }
2327
2328 /**
2329 * Purpose
2330 * =======
2331 *
2332 * DTPMV performs one of the matrix-vector operations
2333 *
2334 * x := A*x, or x := A'*x,
2335 *
2336 * where x is an n element vector and A is an n by n unit, or non-unit,
2337 * upper or lower triangular matrix, supplied in packed form.
2338 *
2339 * Arguments
2340 * ==========
2341 *
2342 * UPLO - CHARACTER*1.
2343 * On entry, UPLO specifies whether the matrix is an upper or
2344 * lower triangular matrix as follows:
2345 *
2346 * UPLO = 'U' or 'u' A is an upper triangular matrix.
2347 *
2348 * UPLO = 'L' or 'l' A is a lower triangular matrix.
2349 *
2350 * Unchanged on exit.
2351 *
2352 * TRANS - CHARACTER*1.
2353 * On entry, TRANS specifies the operation to be performed as
2354 * follows:
2355 *
2356 * TRANS = 'N' or 'n' x := A*x.
2357 *
2358 * TRANS = 'T' or 't' x := A'*x.
2359 *
2360 * TRANS = 'C' or 'c' x := A'*x.
2361 *
2362 * Unchanged on exit.
2363 *
2364 * DIAG - CHARACTER*1.
2365 * On entry, DIAG specifies whether or not A is unit
2366 * triangular as follows:
2367 *
2368 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
2369 *
2370 * DIAG = 'N' or 'n' A is not assumed to be unit
2371 * triangular.
2372 *
2373 * Unchanged on exit.
2374 *
2375 * N - INTEGER.
2376 * On entry, N specifies the order of the matrix A.
2377 * N must be at least zero.
2378 * Unchanged on exit.
2379 *
2380 * AP - DOUBLE PRECISION array of DIMENSION at least
2381 * ( ( n*( n + 1 ) )/2 ).
2382 * Before entry with UPLO = 'U' or 'u', the array AP must
2383 * contain the upper triangular matrix packed sequentially,
2384 * column by column, so that AP( 1 ) contains a( 1, 1 ),
2385 * AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
2386 * respectively, and so on.
2387 * Before entry with UPLO = 'L' or 'l', the array AP must
2388 * contain the lower triangular matrix packed sequentially,
2389 * column by column, so that AP( 1 ) contains a( 1, 1 ),
2390 * AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
2391 * respectively, and so on.
2392 * Note that when DIAG = 'U' or 'u', the diagonal elements of
2393 * A are not referenced, but are assumed to be unity.
2394 * Unchanged on exit.
2395 *
2396 * X - DOUBLE PRECISION array of dimension at least
2397 * ( 1 + ( n - 1 )*abs( INCX ) ).
2398 * Before entry, the incremented array X must contain the n
2399 * element vector x. On exit, X is overwritten with the
2400 * tranformed vector x.
2401 *
2402 * INCX - INTEGER.
2403 * On entry, INCX specifies the increment for the elements of
2404 * X. INCX must not be zero.
2405 * Unchanged on exit.
2406 *
2407 *
2408 * Level 2 Blas routine.
2409 *
2410 * -- Written on 22-October-1986.
2411 * Jack Dongarra, Argonne National Lab.
2412 * Jeremy Du Croz, Nag Central Office.
2413 * Sven Hammarling, Nag Central Office.
2414 * Richard Hanson, Sandia National Labs.
2415 *
2416 *
2417 * .. Parameters ..
2418 **/
C_DTPMV(char uplo,char trans,char diag,int n,double * ap,double * x,int incx)2419 void C_DTPMV(char uplo, char trans, char diag, int n, double *ap, double *x,
2420 int incx)
2421 {
2422 if (n == 0)
2423 return;
2424 if (uplo == 'U' || uplo == 'u')
2425 uplo = 'L';
2426 else if (uplo == 'L' || uplo == 'l')
2427 uplo = 'U';
2428 else
2429 throw std::invalid_argument("C_DTPMV uplo argument is invalid.");
2430 if (trans == 'N' || trans == 'n')
2431 trans = 'T';
2432 else if (trans == 'T' || trans == 't')
2433 trans = 'N';
2434 else
2435 throw std::invalid_argument("C_DTPMV trans argument is invalid.");
2436 ::F_DTPMV(&uplo, &trans, &diag, &n, ap, x, &incx);
2437 }
2438
2439 /**
2440 * Purpose
2441 * =======
2442 *
2443 * DTPSV solves one of the systems of equations
2444 *
2445 * A*x = b, or A'*x = b,
2446 *
2447 * where b and x are n element vectors and A is an n by n unit, or
2448 * non-unit, upper or lower triangular matrix, supplied in packed form.
2449 *
2450 * No test for singularity or near-singularity is included in this
2451 * routine. Such tests must be performed before calling this routine.
2452 *
2453 * Arguments
2454 * ==========
2455 *
2456 * UPLO - CHARACTER*1.
2457 * On entry, UPLO specifies whether the matrix is an upper or
2458 * lower triangular matrix as follows:
2459 *
2460 * UPLO = 'U' or 'u' A is an upper triangular matrix.
2461 *
2462 * UPLO = 'L' or 'l' A is a lower triangular matrix.
2463 *
2464 * Unchanged on exit.
2465 *
2466 * TRANS - CHARACTER*1.
2467 * On entry, TRANS specifies the equations to be solved as
2468 * follows:
2469 *
2470 * TRANS = 'N' or 'n' A*x = b.
2471 *
2472 * TRANS = 'T' or 't' A'*x = b.
2473 *
2474 * TRANS = 'C' or 'c' A'*x = b.
2475 *
2476 * Unchanged on exit.
2477 *
2478 * DIAG - CHARACTER*1.
2479 * On entry, DIAG specifies whether or not A is unit
2480 * triangular as follows:
2481 *
2482 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
2483 *
2484 * DIAG = 'N' or 'n' A is not assumed to be unit
2485 * triangular.
2486 *
2487 * Unchanged on exit.
2488 *
2489 * N - INTEGER.
2490 * On entry, N specifies the order of the matrix A.
2491 * N must be at least zero.
2492 * Unchanged on exit.
2493 *
2494 * AP - DOUBLE PRECISION array of DIMENSION at least
2495 * ( ( n*( n + 1 ) )/2 ).
2496 * Before entry with UPLO = 'U' or 'u', the array AP must
2497 * contain the upper triangular matrix packed sequentially,
2498 * column by column, so that AP( 1 ) contains a( 1, 1 ),
2499 * AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
2500 * respectively, and so on.
2501 * Before entry with UPLO = 'L' or 'l', the array AP must
2502 * contain the lower triangular matrix packed sequentially,
2503 * column by column, so that AP( 1 ) contains a( 1, 1 ),
2504 * AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
2505 * respectively, and so on.
2506 * Note that when DIAG = 'U' or 'u', the diagonal elements of
2507 * A are not referenced, but are assumed to be unity.
2508 * Unchanged on exit.
2509 *
2510 * X - DOUBLE PRECISION array of dimension at least
2511 * ( 1 + ( n - 1 )*abs( INCX ) ).
2512 * Before entry, the incremented array X must contain the n
2513 * element right-hand side vector b. On exit, X is overwritten
2514 * with the solution vector x.
2515 *
2516 * INCX - INTEGER.
2517 * On entry, INCX specifies the increment for the elements of
2518 * X. INCX must not be zero.
2519 * Unchanged on exit.
2520 *
2521 *
2522 * Level 2 Blas routine.
2523 *
2524 * -- Written on 22-October-1986.
2525 * Jack Dongarra, Argonne National Lab.
2526 * Jeremy Du Croz, Nag Central Office.
2527 * Sven Hammarling, Nag Central Office.
2528 * Richard Hanson, Sandia National Labs.
2529 *
2530 *
2531 * .. Parameters ..
2532 **/
C_DTPSV(char uplo,char trans,char diag,int n,double * ap,double * x,int incx)2533 void C_DTPSV(char uplo, char trans, char diag, int n, double *ap, double *x,
2534 int incx)
2535 {
2536 if (n == 0)
2537 return;
2538 if (uplo == 'U' || uplo == 'u')
2539 uplo = 'L';
2540 else if (uplo == 'L' || uplo == 'l')
2541 uplo = 'U';
2542 else
2543 throw std::invalid_argument("C_DTPSV uplo argument is invalid.");
2544 if (trans == 'N' || trans == 'n')
2545 trans = 'T';
2546 else if (trans == 'T' || trans == 't')
2547 trans = 'N';
2548 else
2549 throw std::invalid_argument("C_DTPSV trans argument is invalid.");
2550 ::F_DTPSV(&uplo, &trans, &diag, &n, ap, x, &incx);
2551 }
2552
2553 /**
2554 * Purpose
2555 * =======
2556 *
2557 * DTRMM performs one of the matrix-matrix operations
2558 *
2559 * B := alpha*op( A )*B, or B := alpha*B*op( A ),
2560 *
2561 * where alpha is a scalar, B is an m by n matrix, A is a unit, or
2562 * non-unit, upper or lower triangular matrix and op( A ) is one of
2563 *
2564 * op( A ) = A or op( A ) = A'.
2565 *
2566 * Arguments
2567 * ==========
2568 *
2569 * SIDE - CHARACTER*1.
2570 * On entry, SIDE specifies whether op( A ) multiplies B from
2571 * the left or right as follows:
2572 *
2573 * SIDE = 'L' or 'l' B := alpha*op( A )*B.
2574 *
2575 * SIDE = 'R' or 'r' B := alpha*B*op( A ).
2576 *
2577 * Unchanged on exit.
2578 *
2579 * UPLO - CHARACTER*1.
2580 * On entry, UPLO specifies whether the matrix A is an upper or
2581 * lower triangular matrix as follows:
2582 *
2583 * UPLO = 'U' or 'u' A is an upper triangular matrix.
2584 *
2585 * UPLO = 'L' or 'l' A is a lower triangular matrix.
2586 *
2587 * Unchanged on exit.
2588 *
2589 * TRANSA - CHARACTER*1.
2590 * On entry, TRANSA specifies the form of op( A ) to be used in
2591 * the matrix multiplication as follows:
2592 *
2593 * TRANSA = 'N' or 'n' op( A ) = A.
2594 *
2595 * TRANSA = 'T' or 't' op( A ) = A'.
2596 *
2597 * TRANSA = 'C' or 'c' op( A ) = A'.
2598 *
2599 * Unchanged on exit.
2600 *
2601 * DIAG - CHARACTER*1.
2602 * On entry, DIAG specifies whether or not A is unit triangular
2603 * as follows:
2604 *
2605 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
2606 *
2607 * DIAG = 'N' or 'n' A is not assumed to be unit
2608 * triangular.
2609 *
2610 * Unchanged on exit.
2611 *
2612 * M - INTEGER.
2613 * On entry, M specifies the number of rows of B. M must be at
2614 * least zero.
2615 * Unchanged on exit.
2616 *
2617 * N - INTEGER.
2618 * On entry, N specifies the number of columns of B. N must be
2619 * at least zero.
2620 * Unchanged on exit.
2621 *
2622 * ALPHA - DOUBLE PRECISION.
2623 * On entry, ALPHA specifies the scalar alpha. When alpha is
2624 * zero then A is not referenced and B need not be set before
2625 * entry.
2626 * Unchanged on exit.
2627 *
2628 * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
2629 * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
2630 * Before entry with UPLO = 'U' or 'u', the leading k by k
2631 * upper triangular part of the array A must contain the upper
2632 * triangular matrix and the strictly lower triangular part of
2633 * A is not referenced.
2634 * Before entry with UPLO = 'L' or 'l', the leading k by k
2635 * lower triangular part of the array A must contain the lower
2636 * triangular matrix and the strictly upper triangular part of
2637 * A is not referenced.
2638 * Note that when DIAG = 'U' or 'u', the diagonal elements of
2639 * A are not referenced either, but are assumed to be unity.
2640 * Unchanged on exit.
2641 *
2642 * LDA - INTEGER.
2643 * On entry, LDA specifies the first dimension of A as declared
2644 * in the calling (sub) program. When SIDE = 'L' or 'l' then
2645 * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
2646 * then LDA must be at least max( 1, n ).
2647 * Unchanged on exit.
2648 *
2649 * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
2650 * Before entry, the leading m by n part of the array B must
2651 * contain the matrix B, and on exit is overwritten by the
2652 * transformed matrix.
2653 *
2654 * LDB - INTEGER.
2655 * On entry, LDB specifies the first dimension of B as declared
2656 * in the calling (sub) program. LDB must be at least
2657 * max( 1, m ).
2658 * Unchanged on exit.
2659 *
2660 *
2661 * Level 3 Blas routine.
2662 *
2663 * -- Written on 8-February-1989.
2664 * Jack Dongarra, Argonne National Laboratory.
2665 * Iain Duff, AERE Harwell.
2666 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
2667 * Sven Hammarling, Numerical Algorithms Group Ltd.
2668 *
2669 *
2670 * .. External Functions ..
2671 **/
C_DTRMM(char side,char uplo,char transa,char diag,int m,int n,double alpha,double * a,int lda,double * b,int ldb)2672 void C_DTRMM(char side, char uplo, char transa, char diag, int m, int n,
2673 double alpha, double *a, int lda, double *b, int ldb)
2674 {
2675 if (m == 0 || n == 0)
2676 return;
2677 if (uplo == 'U' || uplo == 'u')
2678 uplo = 'L';
2679 else if (uplo == 'L' || uplo == 'l')
2680 uplo = 'U';
2681 else
2682 throw std::invalid_argument("C_DTRMM uplo argument is invalid.");
2683 if (side == 'L' || side == 'L')
2684 side = 'R';
2685 else if (side == 'R' || side == 'r')
2686 side = 'L';
2687 else
2688 throw std::invalid_argument("C_DTRMM side argument is invalid.");
2689 ::F_DTRMM(&side, &uplo, &transa, &diag, &n, &m, &alpha, a, &lda, b, &ldb);
2690 }
2691
2692 /**
2693 * Purpose
2694 * =======
2695 *
2696 * DTRMV performs one of the matrix-vector operations
2697 *
2698 * x := A*x, or x := A'*x,
2699 *
2700 * where x is an n element vector and A is an n by n unit, or non-unit,
2701 * upper or lower triangular matrix.
2702 *
2703 * Arguments
2704 * ==========
2705 *
2706 * UPLO - CHARACTER*1.
2707 * On entry, UPLO specifies whether the matrix is an upper or
2708 * lower triangular matrix as follows:
2709 *
2710 * UPLO = 'U' or 'u' A is an upper triangular matrix.
2711 *
2712 * UPLO = 'L' or 'l' A is a lower triangular matrix.
2713 *
2714 * Unchanged on exit.
2715 *
2716 * TRANS - CHARACTER*1.
2717 * On entry, TRANS specifies the operation to be performed as
2718 * follows:
2719 *
2720 * TRANS = 'N' or 'n' x := A*x.
2721 *
2722 * TRANS = 'T' or 't' x := A'*x.
2723 *
2724 * TRANS = 'C' or 'c' x := A'*x.
2725 *
2726 * Unchanged on exit.
2727 *
2728 * DIAG - CHARACTER*1.
2729 * On entry, DIAG specifies whether or not A is unit
2730 * triangular as follows:
2731 *
2732 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
2733 *
2734 * DIAG = 'N' or 'n' A is not assumed to be unit
2735 * triangular.
2736 *
2737 * Unchanged on exit.
2738 *
2739 * N - INTEGER.
2740 * On entry, N specifies the order of the matrix A.
2741 * N must be at least zero.
2742 * Unchanged on exit.
2743 *
2744 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
2745 * Before entry with UPLO = 'U' or 'u', the leading n by n
2746 * upper triangular part of the array A must contain the upper
2747 * triangular matrix and the strictly lower triangular part of
2748 * A is not referenced.
2749 * Before entry with UPLO = 'L' or 'l', the leading n by n
2750 * lower triangular part of the array A must contain the lower
2751 * triangular matrix and the strictly upper triangular part of
2752 * A is not referenced.
2753 * Note that when DIAG = 'U' or 'u', the diagonal elements of
2754 * A are not referenced either, but are assumed to be unity.
2755 * Unchanged on exit.
2756 *
2757 * LDA - INTEGER.
2758 * On entry, LDA specifies the first dimension of A as declared
2759 * in the calling (sub) program. LDA must be at least
2760 * max( 1, n ).
2761 * Unchanged on exit.
2762 *
2763 * X - DOUBLE PRECISION array of dimension at least
2764 * ( 1 + ( n - 1 )*abs( INCX ) ).
2765 * Before entry, the incremented array X must contain the n
2766 * element vector x. On exit, X is overwritten with the
2767 * tranformed vector x.
2768 *
2769 * INCX - INTEGER.
2770 * On entry, INCX specifies the increment for the elements of
2771 * X. INCX must not be zero.
2772 * Unchanged on exit.
2773 *
2774 *
2775 * Level 2 Blas routine.
2776 *
2777 * -- Written on 22-October-1986.
2778 * Jack Dongarra, Argonne National Lab.
2779 * Jeremy Du Croz, Nag Central Office.
2780 * Sven Hammarling, Nag Central Office.
2781 * Richard Hanson, Sandia National Labs.
2782 *
2783 *
2784 * .. Parameters ..
2785 **/
C_DTRMV(char uplo,char trans,char diag,int n,double * a,int lda,double * x,int incx)2786 void C_DTRMV(char uplo, char trans, char diag, int n, double *a, int lda,
2787 double *x, int incx)
2788 {
2789 if (n == 0)
2790 return;
2791 if (uplo == 'U' || uplo == 'u')
2792 uplo = 'L';
2793 else if (uplo == 'L' || uplo == 'l')
2794 uplo = 'U';
2795 else
2796 throw std::invalid_argument("C_DTRMV uplo argument is invalid.");
2797 if (trans == 'N' || trans == 'n')
2798 trans = 'T';
2799 else if (trans == 'T' || trans == 't')
2800 trans = 'N';
2801 else
2802 throw std::invalid_argument("C_DTRMV trans argument is invalid.");
2803 ::F_DTRMV(&uplo, &trans, &diag, &n, a, &lda, x, &incx);
2804 }
2805
2806 /**
2807 * Purpose
2808 * =======
2809 *
2810 * DTRSM solves one of the matrix equations
2811 *
2812 * op( A )*X = alpha*B, or X*op( A ) = alpha*B,
2813 *
2814 * where alpha is a scalar, X and B are m by n matrices, A is a unit, or
2815 * non-unit, upper or lower triangular matrix and op( A ) is one of
2816 *
2817 * op( A ) = A or op( A ) = A'.
2818 *
2819 * The matrix X is overwritten on B.
2820 *
2821 * Arguments
2822 * ==========
2823 *
2824 * SIDE - CHARACTER*1.
2825 * On entry, SIDE specifies whether op( A ) appears on the left
2826 * or right of X as follows:
2827 *
2828 * SIDE = 'L' or 'l' op( A )*X = alpha*B.
2829 *
2830 * SIDE = 'R' or 'r' X*op( A ) = alpha*B.
2831 *
2832 * Unchanged on exit.
2833 *
2834 * UPLO - CHARACTER*1.
2835 * On entry, UPLO specifies whether the matrix A is an upper or
2836 * lower triangular matrix as follows:
2837 *
2838 * UPLO = 'U' or 'u' A is an upper triangular matrix.
2839 *
2840 * UPLO = 'L' or 'l' A is a lower triangular matrix.
2841 *
2842 * Unchanged on exit.
2843 *
2844 * TRANSA - CHARACTER*1.
2845 * On entry, TRANSA specifies the form of op( A ) to be used in
2846 * the matrix multiplication as follows:
2847 *
2848 * TRANSA = 'N' or 'n' op( A ) = A.
2849 *
2850 * TRANSA = 'T' or 't' op( A ) = A'.
2851 *
2852 * TRANSA = 'C' or 'c' op( A ) = A'.
2853 *
2854 * Unchanged on exit.
2855 *
2856 * DIAG - CHARACTER*1.
2857 * On entry, DIAG specifies whether or not A is unit triangular
2858 * as follows:
2859 *
2860 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
2861 *
2862 * DIAG = 'N' or 'n' A is not assumed to be unit
2863 * triangular.
2864 *
2865 * Unchanged on exit.
2866 *
2867 * M - INTEGER.
2868 * On entry, M specifies the number of rows of B. M must be at
2869 * least zero.
2870 * Unchanged on exit.
2871 *
2872 * N - INTEGER.
2873 * On entry, N specifies the number of columns of B. N must be
2874 * at least zero.
2875 * Unchanged on exit.
2876 *
2877 * ALPHA - DOUBLE PRECISION.
2878 * On entry, ALPHA specifies the scalar alpha. When alpha is
2879 * zero then A is not referenced and B need not be set before
2880 * entry.
2881 * Unchanged on exit.
2882 *
2883 * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
2884 * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
2885 * Before entry with UPLO = 'U' or 'u', the leading k by k
2886 * upper triangular part of the array A must contain the upper
2887 * triangular matrix and the strictly lower triangular part of
2888 * A is not referenced.
2889 * Before entry with UPLO = 'L' or 'l', the leading k by k
2890 * lower triangular part of the array A must contain the lower
2891 * triangular matrix and the strictly upper triangular part of
2892 * A is not referenced.
2893 * Note that when DIAG = 'U' or 'u', the diagonal elements of
2894 * A are not referenced either, but are assumed to be unity.
2895 * Unchanged on exit.
2896 *
2897 * LDA - INTEGER.
2898 * On entry, LDA specifies the first dimension of A as declared
2899 * in the calling (sub) program. When SIDE = 'L' or 'l' then
2900 * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
2901 * then LDA must be at least max( 1, n ).
2902 * Unchanged on exit.
2903 *
2904 * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
2905 * Before entry, the leading m by n part of the array B must
2906 * contain the right-hand side matrix B, and on exit is
2907 * overwritten by the solution matrix X.
2908 *
2909 * LDB - INTEGER.
2910 * On entry, LDB specifies the first dimension of B as declared
2911 * in the calling (sub) program. LDB must be at least
2912 * max( 1, m ).
2913 * Unchanged on exit.
2914 *
2915 *
2916 * Level 3 Blas routine.
2917 *
2918 *
2919 * -- Written on 8-February-1989.
2920 * Jack Dongarra, Argonne National Laboratory.
2921 * Iain Duff, AERE Harwell.
2922 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
2923 * Sven Hammarling, Numerical Algorithms Group Ltd.
2924 *
2925 *
2926 * .. External Functions ..
2927 **/
C_DTRSM(char side,char uplo,char transa,char diag,int m,int n,double alpha,double * a,int lda,double * b,int ldb)2928 void C_DTRSM(char side, char uplo, char transa, char diag, int m, int n,
2929 double alpha, double *a, int lda, double *b, int ldb)
2930 {
2931 if (m == 0 || n == 0)
2932 return;
2933 if (uplo == 'U' || uplo == 'u')
2934 uplo = 'L';
2935 else if (uplo == 'L' || uplo == 'l')
2936 uplo = 'U';
2937 else
2938 throw std::invalid_argument("C_DTRSM uplo argument is invalid.");
2939 if (side == 'L' || side == 'L')
2940 side = 'R';
2941 else if (side == 'R' || side == 'r')
2942 side = 'L';
2943 else
2944 throw std::invalid_argument("C_DTRSM side argument is invalid.");
2945 ::F_DTRSM(&side, &uplo, &transa, &diag, &n, &m, &alpha, a, &lda, b, &ldb);
2946 }
2947
2948 /**
2949 * Purpose
2950 * =======
2951 *
2952 * DTRSV solves one of the systems of equations
2953 *
2954 * A*x = b, or A'*x = b,
2955 *
2956 * where b and x are n element vectors and A is an n by n unit, or
2957 * non-unit, upper or lower triangular matrix.
2958 *
2959 * No test for singularity or near-singularity is included in this
2960 * routine. Such tests must be performed before calling this routine.
2961 *
2962 * Arguments
2963 * ==========
2964 *
2965 * UPLO - CHARACTER*1.
2966 * On entry, UPLO specifies whether the matrix is an upper or
2967 * lower triangular matrix as follows:
2968 *
2969 * UPLO = 'U' or 'u' A is an upper triangular matrix.
2970 *
2971 * UPLO = 'L' or 'l' A is a lower triangular matrix.
2972 *
2973 * Unchanged on exit.
2974 *
2975 * TRANS - CHARACTER*1.
2976 * On entry, TRANS specifies the equations to be solved as
2977 * follows:
2978 *
2979 * TRANS = 'N' or 'n' A*x = b.
2980 *
2981 * TRANS = 'T' or 't' A'*x = b.
2982 *
2983 * TRANS = 'C' or 'c' A'*x = b.
2984 *
2985 * Unchanged on exit.
2986 *
2987 * DIAG - CHARACTER*1.
2988 * On entry, DIAG specifies whether or not A is unit
2989 * triangular as follows:
2990 *
2991 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
2992 *
2993 * DIAG = 'N' or 'n' A is not assumed to be unit
2994 * triangular.
2995 *
2996 * Unchanged on exit.
2997 *
2998 * N - INTEGER.
2999 * On entry, N specifies the order of the matrix A.
3000 * N must be at least zero.
3001 * Unchanged on exit.
3002 *
3003 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
3004 * Before entry with UPLO = 'U' or 'u', the leading n by n
3005 * upper triangular part of the array A must contain the upper
3006 * triangular matrix and the strictly lower triangular part of
3007 * A is not referenced.
3008 * Before entry with UPLO = 'L' or 'l', the leading n by n
3009 * lower triangular part of the array A must contain the lower
3010 * triangular matrix and the strictly upper triangular part of
3011 * A is not referenced.
3012 * Note that when DIAG = 'U' or 'u', the diagonal elements of
3013 * A are not referenced either, but are assumed to be unity.
3014 * Unchanged on exit.
3015 *
3016 * LDA - INTEGER.
3017 * On entry, LDA specifies the first dimension of A as declared
3018 * in the calling (sub) program. LDA must be at least
3019 * max( 1, n ).
3020 * Unchanged on exit.
3021 *
3022 * X - DOUBLE PRECISION array of dimension at least
3023 * ( 1 + ( n - 1 )*abs( INCX ) ).
3024 * Before entry, the incremented array X must contain the n
3025 * element right-hand side vector b. On exit, X is overwritten
3026 * with the solution vector x.
3027 *
3028 * INCX - INTEGER.
3029 * On entry, INCX specifies the increment for the elements of
3030 * X. INCX must not be zero.
3031 * Unchanged on exit.
3032 *
3033 *
3034 * Level 2 Blas routine.
3035 *
3036 * -- Written on 22-October-1986.
3037 * Jack Dongarra, Argonne National Lab.
3038 * Jeremy Du Croz, Nag Central Office.
3039 * Sven Hammarling, Nag Central Office.
3040 * Richard Hanson, Sandia National Labs.
3041 *
3042 *
3043 * .. Parameters ..
3044 **/
C_DTRSV(char uplo,char trans,char diag,int n,double * a,int lda,double * x,int incx)3045 void C_DTRSV(char uplo, char trans, char diag, int n, double *a, int lda,
3046 double *x, int incx)
3047 {
3048 if (n == 0)
3049 return;
3050 if (uplo == 'U' || uplo == 'u')
3051 uplo = 'L';
3052 else if (uplo == 'L' || uplo == 'l')
3053 uplo = 'U';
3054 else
3055 throw std::invalid_argument("C_DTRSV uplo argument is invalid.");
3056 if (trans == 'N' || trans == 'n')
3057 trans = 'T';
3058 else if (trans == 'T' || trans == 't')
3059 trans = 'N';
3060 else
3061 throw std::invalid_argument("C_DTRSV trans argument is invalid.");
3062 ::F_DTRSV(&uplo, &trans, &diag, &n, a, &lda, x, &incx);
3063 }
3064 }
3065