ergo
template_blas_ger.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_BLAS_GER_HEADER
36 #define TEMPLATE_BLAS_GER_HEADER
37 
38 
39 template<class Treal>
40 int template_blas_ger(const integer *m, const integer *n, const Treal *alpha,
41  const Treal *x, const integer *incx, const Treal *y, const integer *incy,
42  Treal *a, const integer *lda)
43 {
44  /* System generated locals */
45  integer a_dim1, a_offset, i__1, i__2;
46  /* Local variables */
47  integer info;
48  Treal temp;
49  integer i__, j, ix, jy, kx;
50 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
51 /* Purpose
52  =======
53  DGER performs the rank 1 operation
54  A := alpha*x*y' + A,
55  where alpha is a scalar, x is an m element vector, y is an n element
56  vector and A is an m by n matrix.
57  Parameters
58  ==========
59  M - INTEGER.
60  On entry, M specifies the number of rows of the matrix A.
61  M must be at least zero.
62  Unchanged on exit.
63  N - INTEGER.
64  On entry, N specifies the number of columns of the matrix A.
65  N must be at least zero.
66  Unchanged on exit.
67  ALPHA - DOUBLE PRECISION.
68  On entry, ALPHA specifies the scalar alpha.
69  Unchanged on exit.
70  X - DOUBLE PRECISION array of dimension at least
71  ( 1 + ( m - 1 )*abs( INCX ) ).
72  Before entry, the incremented array X must contain the m
73  element vector x.
74  Unchanged on exit.
75  INCX - INTEGER.
76  On entry, INCX specifies the increment for the elements of
77  X. INCX must not be zero.
78  Unchanged on exit.
79  Y - DOUBLE PRECISION array of dimension at least
80  ( 1 + ( n - 1 )*abs( INCY ) ).
81  Before entry, the incremented array Y must contain the n
82  element vector y.
83  Unchanged on exit.
84  INCY - INTEGER.
85  On entry, INCY specifies the increment for the elements of
86  Y. INCY must not be zero.
87  Unchanged on exit.
88  A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
89  Before entry, the leading m by n part of the array A must
90  contain the matrix of coefficients. On exit, A is
91  overwritten by the updated matrix.
92  LDA - INTEGER.
93  On entry, LDA specifies the first dimension of A as declared
94  in the calling (sub) program. LDA must be at least
95  max( 1, m ).
96  Unchanged on exit.
97  Level 2 Blas routine.
98  -- Written on 22-October-1986.
99  Jack Dongarra, Argonne National Lab.
100  Jeremy Du Croz, Nag Central Office.
101  Sven Hammarling, Nag Central Office.
102  Richard Hanson, Sandia National Labs.
103  Test the input parameters.
104  Parameter adjustments */
105  --x;
106  --y;
107  a_dim1 = *lda;
108  a_offset = 1 + a_dim1 * 1;
109  a -= a_offset;
110  /* Function Body */
111  info = 0;
112  if (*m < 0) {
113  info = 1;
114  } else if (*n < 0) {
115  info = 2;
116  } else if (*incx == 0) {
117  info = 5;
118  } else if (*incy == 0) {
119  info = 7;
120  } else if (*lda < maxMACRO(1,*m)) {
121  info = 9;
122  }
123  if (info != 0) {
124  template_blas_erbla("GER ", &info);
125  return 0;
126  }
127 /* Quick return if possible. */
128  if (*m == 0 || *n == 0 || *alpha == 0.) {
129  return 0;
130  }
131 /* Start the operations. In this version the elements of A are
132  accessed sequentially with one pass through A. */
133  if (*incy > 0) {
134  jy = 1;
135  } else {
136  jy = 1 - (*n - 1) * *incy;
137  }
138  if (*incx == 1) {
139  i__1 = *n;
140  for (j = 1; j <= i__1; ++j) {
141  if (y[jy] != 0.) {
142  temp = *alpha * y[jy];
143  i__2 = *m;
144  for (i__ = 1; i__ <= i__2; ++i__) {
145  a_ref(i__, j) = a_ref(i__, j) + x[i__] * temp;
146 /* L10: */
147  }
148  }
149  jy += *incy;
150 /* L20: */
151  }
152  } else {
153  if (*incx > 0) {
154  kx = 1;
155  } else {
156  kx = 1 - (*m - 1) * *incx;
157  }
158  i__1 = *n;
159  for (j = 1; j <= i__1; ++j) {
160  if (y[jy] != 0.) {
161  temp = *alpha * y[jy];
162  ix = kx;
163  i__2 = *m;
164  for (i__ = 1; i__ <= i__2; ++i__) {
165  a_ref(i__, j) = a_ref(i__, j) + x[ix] * temp;
166  ix += *incx;
167 /* L30: */
168  }
169  }
170  jy += *incy;
171 /* L40: */
172  }
173  }
174  return 0;
175 /* End of DGER . */
176 } /* dger_ */
177 #undef a_ref
178 
179 #endif