2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 #include "gmx_lapack.h"
38 #include "lapack_limits.h"
42 F77_FUNC(slatrd,SLATRD)(const char * uplo,
54 float one,zero,minusone,alpha;
55 const char ch=toupper(*uplo);
65 for(i=*n;i>=(*n-*nb+1);i--) {
72 F77_FUNC(sgemv,SGEMV)("N",&i,&ti1,&minusone, &(a[ i*(*lda) + 0]),lda,&(w[iw*(*ldw)+(i-1)]),
73 ldw,&one, &(a[ (i-1)*(*lda) + 0]), &ti2);
75 F77_FUNC(sgemv,SGEMV)("N",&i,&ti1,&minusone, &(w[ iw*(*ldw) + 0]),ldw,&(a[i*(*lda)+(i-1)]),
76 lda,&one, &(a[ (i-1)*(*lda) + 0]), &ti2);
80 /* Generate elementary reflector H(i) to annihilate
87 F77_FUNC(slarfg,SLARFG)(&ti1,&(a[(i-1)*(*lda)+(i-2)]),&(a[(i-1)*(*lda)+0]),&ti2,&(tau[i-2]));
89 e[i-2] = a[(i-1)*(*lda)+(i-2)];
90 a[(i-1)*(*lda)+(i-2)] = 1.0;
92 /* Compute W(1:i-1,i) */
97 F77_FUNC(ssymv,SSYMV)("U",&ti1,&one,a,lda,&(a[(i-1)*(*lda)+0]),&ti2,&zero,
98 &(w[(iw-1)*(*ldw)+0]),&ti2);
104 F77_FUNC(sgemv,SGEMV)("T",&ti1,&ti2,&one,&(w[iw*(*ldw)+0]),ldw,&(a[(i-1)*(*lda)+0]),&ti3,
105 &zero,&(w[(iw-1)*(*ldw)+i]),&ti3);
108 F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone,&(a[i*(*lda)+0]),lda,&(w[(iw-1)*(*ldw)+i]),&ti3,
109 &one,&(w[(iw-1)*(*ldw)+0]),&ti3);
112 F77_FUNC(sgemv,SGEMV)("T",&ti1,&ti2,&one,&(a[i*(*lda)+0]),lda,&(a[(i-1)*(*lda)+0]),&ti3,
113 &zero,&(w[(iw-1)*(*ldw)+i]),&ti3);
116 F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone,&(w[iw*(*ldw)+0]),ldw,&(w[(iw-1)*(*ldw)+i]),&ti3,
117 &one,&(w[(iw-1)*(*ldw)+0]),&ti3);
123 F77_FUNC(sscal,SSCAL)(&ti1,&(tau[i-2]),&(w[(iw-1)*(*ldw)+0]),&ti2);
125 alpha = -0.5*tau[i-2]*F77_FUNC(sdot,SDOT)(&ti1,&(w[(iw-1)*(*ldw)+0]),&ti2,
126 &(a[(i-1)*(*lda)+0]),&ti2);
129 F77_FUNC(saxpy,SAXPY)(&ti1,&alpha,&(a[(i-1)*(*lda)+0]),&ti2,&(w[(iw-1)*(*ldw)+0]),&ti2);
135 for(i=1;i<=*nb;i++) {
141 F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone, &(a[ i-1 ]),lda,&(w[ i-1 ]),
142 ldw,&one, &(a[ (i-1)*(*lda) + (i-1)]), &ti3);
144 F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone, &(w[ i-1 ]),ldw,&(a[ i-1 ]),
145 lda,&one, &(a[ (i-1)*(*lda) + (i-1)]), &ti3);
149 ti2 = (*n < i+2 ) ? *n : (i+2);
152 F77_FUNC(slarfg,SLARFG)(&ti1,&(a[(i-1)*(*lda)+(i)]),&(a[(i-1)*(*lda)+(ti2-1)]),&ti3,&(tau[i-1]));
153 e[i-1] = a[(i-1)*(*lda)+(i)];
154 a[(i-1)*(*lda)+(i)] = 1.0;
158 F77_FUNC(ssymv,SSYMV)("L",&ti1,&one,&(a[i*(*lda)+i]),lda,&(a[(i-1)*(*lda)+i]),&ti2,
159 &zero,&(w[(i-1)*(*ldw)+i]),&ti2);
164 F77_FUNC(sgemv,SGEMV)("T",&ti1,&ti2,&one,&(w[ i ]),ldw,&(a[(i-1)*(*lda)+i]),&ti3,
165 &zero,&(w[(i-1)*(*ldw)+0]),&ti3);
168 F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone,&(a[ i ]),lda,&(w[(i-1)*(*ldw)+0]),&ti3,
169 &one,&(w[(i-1)*(*ldw)+i]),&ti3);
172 F77_FUNC(sgemv,SGEMV)("T",&ti1,&ti2,&one,&(a[ i ]),lda,&(a[(i-1)*(*lda)+i]),&ti3,
173 &zero,&(w[(i-1)*(*ldw)+0]),&ti3);
176 F77_FUNC(sgemv,SGEMV)("N",&ti1,&ti2,&minusone,&(w[ i ]),ldw,&(w[(i-1)*(*ldw)+0]),&ti3,
177 &one,&(w[(i-1)*(*ldw)+i]),&ti3);
179 F77_FUNC(sscal,SSCAL)(&ti1,&(tau[i-1]),&(w[(i-1)*(*ldw)+i]),&ti3);
180 alpha = -0.5*tau[i-1]*F77_FUNC(sdot,SDOT)(&ti1,&(w[(i-1)*(*ldw)+i]),&ti3,
181 &(a[(i-1)*(*lda)+i]),&ti3);
183 F77_FUNC(saxpy,SAXPY)(&ti1,&alpha,&(a[(i-1)*(*lda)+i]),&ti3,&(w[(i-1)*(*ldw)+i]),&ti3);