1
/*! \defgroup QT libqt: The Quantum-Trio Miscellaneous Library */
5
\brief The PSI3 BLAS interface routines
7
Interface to the BLAS routines
14
Additions by TD Crawford and EF Valeev, June 1999.
20
#define F_DAXPY daxpy_
21
#define F_DCOPY dcopy_
22
#define F_DGEMM dgemm_
24
#define F_DSCAL dscal_
25
#define F_DGEMV dgemv_
26
#define F_DSPMV dspmv_
47
#define F_DAXPY DAXPY_
48
#define F_DCOPY DCOPY_
49
#define F_DGEMM DGEMM_
51
#define F_DSCAL DSCAL_
52
#define F_DGEMV DGEMV_
53
#define F_DSPMV DSPMV_
57
extern void F_DAXPY(int *length, double *a, double *x, int *inc_x,
58
double *y, int *inc_y);
59
extern void F_DCOPY(int *length, double *x, int *inc_x,
60
double *y, int *inc_y);
61
extern void F_DGEMM(char *transa, char *transb, int *m, int *n, int *k,
62
double *alpha, double *A, int *lda, double *B, int *ldb,
63
double *beta, double *C, int *ldc);
64
extern void F_DROT(int *ntot,double *x, int *incx,double *y, int *incy,
65
double *cotheta,double *sintheta);
66
extern void F_DSCAL(int *n, double *alpha, double *vec, int *inc);
67
extern void F_DGEMV(char *transa, int *m, int *n, double *alpha, double *A,
68
int *lda, double *X, int *inc_x, double *beta,
69
double *Y, int *inc_y);
70
extern double F_DDOT(int *n, double *x, int *incx, double *y, int *incy);
74
** This function performs y = a * x + y.
75
** Steps every inc_x in x and every inc_y in y (normally both 1).
78
void C_DAXPY(int length, double a, double *x, int inc_x,
81
F_DAXPY(&length, &a, x, &inc_x, y, &inc_y);
86
** This function copies x into y.
87
** Steps every inc_x in x and every inc_y in y (normally both 1).
90
void C_DCOPY(int length, double *x, int inc_x,
93
F_DCOPY(&length, x, &inc_x, y, &inc_y);
99
** This function scales a vector by a real scalar.
102
void C_DSCAL(int n, double alpha, double *vec, int inc)
104
F_DSCAL(&n, &alpha, vec, &inc);
110
** This function calculates plane Givens rotation for vectors
111
** x,y and angle theta. x = x*cos + y*sin, y = -x*sin + y*cos.
113
** \param x: vector x
114
** \param y: vector Y
115
** \param tot: length of x,y
116
** \param incx: increment for x
117
** \param incy: increment for y
120
void C_DROT(int ntot, double *x, int incx, double *y, int incy,
121
double costheta, double sintheta)
124
F_DROT(&ntot,x,&incx,y,&incy,&costheta,&sintheta);
131
** This function calculates C(m,n)=alpha*(opT)A(m,k)*(opT)B(k,n)+ beta*C(m,n)
133
** These arguments mimic their Fortran conterparts; parameters have been
134
** reversed nca, ncb, ncc, A, B, C, to make it correct for C.
136
** \param char transa: On entry, specifies the form of (op)A used in the
137
** matrix multiplication:
138
** If transa = 'N' or 'n', (op)A = A.
139
** If transa = 'T' or 't', (op)A = transp(A).
140
** If transa = 'R' or 'r', (op)A = conjugate(A).
141
** If transa = 'C' or 'c', (op)A = conjug_transp(A).
142
** On exit, transa is unchanged.
144
** \param char transb: On entry, specifies the form of (op)B used in the
145
** matrix multiplication:
146
** If transb = 'N' or 'n', (op)B = B.
147
** If transb = 'T' or 't', (op)B = transp(B).
148
** If transb = 'R' or 'r', (op)B = conjugate(B)
150
** \param int m: On entry, the number of rows of the matrix (op)A and of
151
** the matrix C; m >= 0. On exit, m is unchanged.
153
** \param int n: On entry, the number of columns of the matrix (op)B and
154
** of the matrix C; n >= 0. On exit, n is unchanged.
156
** \param int k: On entry, the number of columns of the matrix (op)A and
157
** the number of rows of the matrix (op)B; k >= 0. On exit,
160
** \param double alpha: On entry, specifies the scalar alpha. On exit, alpha is
163
** \param double* A: On entry, a two-dimensional array A with dimensions ka
164
** by nca. For (op)A = A or conjugate(A), nca >= k and the
165
** leading m by k portion of the array A contains the matrix
166
** A. For (op)A = transp(A) or conjug_transp(A), nca >= m
167
** and the leading k by m part of the array A contains the
168
** matrix A. On exit, a is unchanged.
170
** \param int nca: On entry, the second dimension of array A.
171
** For (op)A = A or conjugate(A), nca >= MAX(1,k).
172
** For (op)A=transp(A) or conjug_transp(A), nca >= MAX(1,m).
173
** On exit, nca is unchanged.
175
** \param double* B: On entry, a two-dimensional array B with dimensions kb
176
** by ncb. For (op)B = B or conjugate(B), kb >= k and the
177
** leading k by n portion of the array contains the matrix
178
** B. For (op)B = transp(B) or conjug_transp(B), ncb >= k and
179
** the leading n by k part of the array contains the matrix
180
** B. On exit, B is unchanged.
182
** \param int ncb: On entry, the second dimension of array B.
183
** For (op)B = B or <conjugate(B), ncb >= MAX(1,n).
184
** For (op)B = transp(B) or conjug_transp(B), ncb >=
185
** MAX(1,k). On exit, ncb is unchanged.
187
** \param double beta: On entry, specifies the scalar beta. On exit, beta is
190
** \param double* C: On entry, a two-dimensional array with the dimension
191
** at least m by ncc. On exit, the leading m by n part of
192
** array C is overwritten by the matrix alpha*(op)A*(op)B +
195
** \param int ncc: On entry, the second dimension of array C;
196
** ncc >=MAX(1,n). On exit, ncc is unchanged.
199
void C_DGEMM(char transa, char transb, int m, int n, int k, double alpha,
200
double *A, int nca, double *B, int ncb, double beta, double *C,
204
/* the only strange thing we need to do is reverse everything
205
since the stride runs differently in C vs. Fortran
208
/* also, do nothing if a dimension is 0 */
209
if (m == 0 || n == 0 || k == 0) return;
211
F_DGEMM(&transb,&transa,&n,&m,&k,&alpha,B,&ncb,A,&nca,&beta,C,&ncc);
217
** This function calculates the matrix-vector product:
219
** Y = alpha * A * X + beta * Y
221
** where X and Y are vectors, A is a matrix, and alpha and beta are
224
** \param char transa: Indicates whether the matrix A should be
225
** transposed ('t') or left alone ('n').
226
** \param int m: The row dimension of A (regardless of transa).
227
** \param int n: The column dimension of A (regardless of transa).
228
** \param double alpha: The scalar alpha.
229
** \param double* A: A pointer to the beginning of the data in A.
230
** \param int nca: The number of columns *actually* in A. This is
231
** useful if one only wishes to multiply the first
232
** n columns of A times X even though A
233
** contains nca columns.
234
** \param double* X: A pointer to the beginning of the data in X.
235
** \param int inc_x: The desired stride for X. Useful for skipping
236
** sections of data to treat only one column of a
237
** complete matrix. Usually 1, though.
238
** \param double beta: The scalar beta.
239
** \param double* Y: A pointer to the beginning of the data in Y.
240
** \param int inc_y: The desired stride for Y.
242
** Interface written by TD Crawford and EF Valeev.
247
void C_DGEMV(char transa, int m, int n, double alpha, double *A,
248
int nca, double *X, int inc_x, double beta, double *Y,
251
if (m == 0 || n == 0) return;
253
if(transa == 'n') transa = 't';
256
F_DGEMV(&transa,&n,&m,&alpha,A,&nca,X,&inc_x,&beta,Y,&inc_y);
263
** This function calculates the matrix-vector product:
265
** Y = alpha * A * X + beta * Y
267
** where X and Y are vectors, A is a matrix, and alpha and beta are
270
** \param char uplo: Indicates whether the matrix A is packed in
271
** upper ('U' or 'u') or lower ('L' or 'l')
272
** triangular form. We reverse what is passed
273
** before sending it on to Fortran because of
274
** the different Fortran/C conventions
275
** \param int n: The order of the matrix A (number of rows/columns)
276
** \param double alpha: The scalar alpha.
277
** \param double* A: A pointer to the beginning of the data in A.
278
** \param double* X: A pointer to the beginning of the data in X.
279
** \param int inc_x: The desired stride for X. Useful for skipping
280
** sections of data to treat only one column of a
281
** complete matrix. Usually 1, though.
282
** \param double beta: The scalar beta.
283
** \param double* Y: A pointer to the beginning of the data in Y.
284
** \param int inc_y: The desired stride for Y.
286
** Interface written by CD Sherrill
291
void C_DSPMV(char uplo, int n, double alpha, double *A,
292
double *X, int inc_x, double beta, double *Y,
297
if (uplo != 'U' && uplo != 'u' && uplo != 'L' && uplo != 'l')
298
fprintf(stderr, "C_DSPMV: called with unrecognized option for uplo!\n");
300
if (uplo == 'U' || uplo == 'u') uplo = 'L';
303
F_DSPMV(&uplo,&n,&alpha,A,X,&inc_x,&beta,Y,&inc_y);
311
** This function returns the dot product of two vectors, X and Y.
313
** \param n: Number of elements in X and Y.
314
** \param X: A pointer to the beginning of the data in X.
315
** Must be of at least length (1+(N-1)*abs(inc_x).
316
** \param inc_x: The desired stride of X. Useful for skipping
317
** around to certain data stored in X.
318
** \param Y: A pointer to the beginning of the data in Y.
319
** \param inc_y: The desired stride for Y.
321
** Interface written by ST Brown.
326
double C_DDOT(int n, double *X, int inc_x, double *Y, int inc_y)
328
if(n == 0) return 0.0;
330
return F_DDOT(&n,X,&inc_x,Y,&inc_y);