~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/lib/libqt/blas_intfc.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
/*! \defgroup QT libqt: The Quantum-Trio Miscellaneous Library */
2
 
 
3
 
/*! \file blas_intfc.c
4
 
    \ingroup (QT)
5
 
    \brief The PSI3 BLAS interface routines
6
 
 
7
 
 Interface to the BLAS routines
8
 
 
9
 
 C. David Sherrill
10
 
 Anna I. Krylov
11
 
 
12
 
 May 1998
13
 
 
14
 
 Additions by TD Crawford and EF Valeev, June 1999.
15
 
*/
16
 
 
17
 
#include <stdio.h>
18
 
 
19
 
#if FC_SYMBOL==2
20
 
#define F_DAXPY daxpy_
21
 
#define F_DCOPY dcopy_
22
 
#define F_DGEMM dgemm_
23
 
#define F_DROT drot_
24
 
#define F_DSCAL dscal_
25
 
#define F_DGEMV dgemv_
26
 
#define F_DSPMV dspmv_
27
 
#define F_DDOT  ddot_
28
 
#elif FC_SYMBOL==1
29
 
#define F_DAXPY daxpy
30
 
#define F_DCOPY dcopy
31
 
#define F_DGEMM dgemm
32
 
#define F_DROT drot
33
 
#define F_DSCAL dscal
34
 
#define F_DGEMV dgemv
35
 
#define F_DSPMV dspmv
36
 
#define F_DDOT  ddot
37
 
#elif FC_SYMBOL==3
38
 
#define F_DAXPY DAXPY
39
 
#define F_DCOPY DCOPY
40
 
#define F_DGEMM DGEMM
41
 
#define F_DROT DROT
42
 
#define F_DSCAL DSCAL
43
 
#define F_DGEMV DGEMV
44
 
#define F_DSPMV DSPMV
45
 
#define F_DDOT  DDOT
46
 
#elif FC_SYMBOL==4
47
 
#define F_DAXPY DAXPY_
48
 
#define F_DCOPY DCOPY_
49
 
#define F_DGEMM DGEMM_
50
 
#define F_DROT DROT_
51
 
#define F_DSCAL DSCAL_
52
 
#define F_DGEMV DGEMV_
53
 
#define F_DSPMV DSPMV_
54
 
#define F_DDOT  DDOT_
55
 
#endif
56
 
 
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);
71
 
 
72
 
/*! 
73
 
** C_DAXPY()
74
 
** This function performs y = a * x + y. 
75
 
** Steps every inc_x in x and every inc_y in y (normally both 1).
76
 
** \ingroup (QT)
77
 
*/
78
 
void C_DAXPY(int length, double a, double *x, int inc_x, 
79
 
             double *y, int inc_y)
80
 
{
81
 
  F_DAXPY(&length, &a, x, &inc_x, y, &inc_y);
82
 
}
83
 
 
84
 
/*!
85
 
** C_DCOPY()
86
 
** This function copies x into y.
87
 
** Steps every inc_x in x and every inc_y in y (normally both 1).
88
 
** \ingroup (QT)
89
 
*/
90
 
void C_DCOPY(int length, double *x, int inc_x, 
91
 
             double *y, int inc_y)
92
 
{
93
 
  F_DCOPY(&length, x, &inc_x, y, &inc_y);
94
 
}
95
 
 
96
 
 
97
 
/*!
98
 
** C_DSCAL()
99
 
** This function scales a vector by a real scalar.
100
 
** \ingroup (QT)
101
 
*/
102
 
void C_DSCAL(int n, double alpha, double *vec, int inc)
103
 
{
104
 
  F_DSCAL(&n, &alpha, vec, &inc);
105
 
}
106
 
 
107
 
 
108
 
/*!
109
 
** C_DROT()
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.
112
 
**
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
118
 
** \ingroup (QT)
119
 
*/
120
 
void C_DROT(int ntot, double *x, int incx, double *y, int incy,
121
 
            double costheta, double sintheta)
122
 
{
123
 
 
124
 
  F_DROT(&ntot,x,&incx,y,&incy,&costheta,&sintheta);
125
 
}
126
 
 
127
 
 
128
 
/*!
129
 
** C_DGEMM()
130
 
**
131
 
** This function calculates C(m,n)=alpha*(opT)A(m,k)*(opT)B(k,n)+ beta*C(m,n)
132
 
**
133
 
** These arguments mimic their Fortran conterparts; parameters have been 
134
 
** reversed nca, ncb, ncc, A, B, C,  to make it correct for C.
135
 
**  
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.
143
 
**
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)
149
 
** 
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.
152
 
**
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.
155
 
**
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,
158
 
**                    k is unchanged.
159
 
** 
160
 
** \param double alpha:  On entry, specifies the scalar alpha. On exit, alpha is
161
 
**                       unchanged.
162
 
** 
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.
169
 
** 
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.
174
 
**
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.
181
 
**
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.
186
 
**
187
 
** \param double beta: On entry, specifies the scalar beta. On exit, beta is
188
 
**                    unchanged.
189
 
**
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 +
193
 
**                    beta*C.
194
 
**
195
 
** \param int ncc:    On entry, the second dimension  of array C; 
196
 
**                    ncc >=MAX(1,n).  On exit, ncc is unchanged.
197
 
** \ingroup (QT)
198
 
*/
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,
201
 
           int ncc)
202
 
{
203
 
 
204
 
  /* the only strange thing we need to do is reverse everything
205
 
     since the stride runs differently in C vs. Fortran
206
 
   */
207
 
  
208
 
  /* also, do nothing if a dimension is 0 */
209
 
  if (m == 0 || n == 0 || k == 0) return;
210
 
 
211
 
  F_DGEMM(&transb,&transa,&n,&m,&k,&alpha,B,&ncb,A,&nca,&beta,C,&ncc);
212
 
 
213
 
}
214
 
 
215
 
/*!
216
 
**  C_DGEMV()
217
 
**  This function calculates the matrix-vector product:
218
 
**
219
 
**  Y = alpha * A * X + beta * Y
220
 
**
221
 
** where X and Y are vectors, A is a matrix, and alpha and beta are
222
 
** constants. 
223
 
**
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.
241
 
** 
242
 
** Interface written by TD Crawford and EF Valeev.
243
 
** June 1999.
244
 
** \ingroup (QT)
245
 
*/
246
 
 
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,
249
 
             int inc_y)
250
 
{
251
 
  if (m == 0 || n == 0) return;
252
 
 
253
 
  if(transa == 'n') transa = 't';
254
 
  else transa = 'n';
255
 
 
256
 
  F_DGEMV(&transa,&n,&m,&alpha,A,&nca,X,&inc_x,&beta,Y,&inc_y);
257
 
 
258
 
}
259
 
 
260
 
 
261
 
/*!
262
 
**  C_DSPMV()
263
 
**  This function calculates the matrix-vector product:
264
 
**
265
 
**  Y = alpha * A * X + beta * Y
266
 
**
267
 
** where X and Y are vectors, A is a matrix, and alpha and beta are
268
 
** constants. 
269
 
**
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.
285
 
** 
286
 
** Interface written by CD Sherrill
287
 
** July 2003
288
 
** \ingroup (QT)
289
 
*/
290
 
 
291
 
void C_DSPMV(char uplo, int n, double alpha, double *A, 
292
 
             double *X, int inc_x, double beta, double *Y,
293
 
             int inc_y)
294
 
{
295
 
  if (n == 0) return;
296
 
 
297
 
  if (uplo != 'U' && uplo != 'u' && uplo != 'L' && uplo != 'l')
298
 
    fprintf(stderr, "C_DSPMV: called with unrecognized option for uplo!\n");
299
 
 
300
 
  if (uplo == 'U' || uplo == 'u') uplo = 'L';
301
 
  else uplo = 'U';
302
 
 
303
 
  F_DSPMV(&uplo,&n,&alpha,A,X,&inc_x,&beta,Y,&inc_y);
304
 
 
305
 
}
306
 
 
307
 
 
308
 
 
309
 
/*!
310
 
** C_DDOT()
311
 
** This function returns the dot product of two vectors, X and Y.
312
 
**
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.
320
 
**
321
 
** Interface written by ST Brown.
322
 
** July 2000
323
 
** \ingroup (QT)
324
 
*/
325
 
 
326
 
double C_DDOT(int n, double *X, int inc_x, double *Y, int inc_y)
327
 
{
328
 
   if(n == 0) return 0.0;
329
 
 
330
 
   return F_DDOT(&n,X,&inc_x,Y,&inc_y);
331
 
}
332