2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012, 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 <types/simple.h>
42 F77_FUNC(dtrmm,DTRMM)(const char *side,
54 int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
60 double alpha = *alpha__;
70 a_offset = 1 + a_dim1;
73 b_offset = 1 + b_dim1;
77 lside = (*side=='L' || *side=='l');
83 nounit = (*diag=='N' || *diag=='n');
84 upper = (*uplo=='U' || *uplo=='u');
91 if (fabs(alpha)<GMX_DOUBLE_MIN) {
93 for (j = 1; j <= i__1; ++j) {
95 for (i__ = 1; i__ <= i__2; ++i__) {
96 b[i__ + j * b_dim1] = 0.;
102 if (*transa=='N' || *transa=='n') {
105 for (j = 1; j <= i__1; ++j) {
107 for (k = 1; k <= i__2; ++k) {
108 if (fabs(b[k + j * b_dim1])>GMX_DOUBLE_MIN) {
109 temp = alpha * b[k + j * b_dim1];
111 for (i__ = 1; i__ <= i__3; ++i__) {
112 b[i__ + j * b_dim1] += temp * a[i__ + k * a_dim1];
115 temp *= a[k + k * a_dim1];
117 b[k + j * b_dim1] = temp;
123 for (j = 1; j <= i__1; ++j) {
124 for (k = m; k >= 1; --k) {
125 if (fabs(b[k + j * b_dim1])>GMX_DOUBLE_MIN) {
126 temp = alpha * b[k + j * b_dim1];
127 b[k + j * b_dim1] = temp;
129 b[k + j * b_dim1] *= a[k + k * a_dim1];
132 for (i__ = k + 1; i__ <= i__2; ++i__) {
133 b[i__ + j * b_dim1] += temp * a[i__ + k *
144 for (j = 1; j <= i__1; ++j) {
145 for (i__ = m; i__ >= 1; --i__) {
146 temp = b[i__ + j * b_dim1];
148 temp *= a[i__ + i__ * a_dim1];
151 for (k = 1; k <= i__2; ++k) {
152 temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
154 b[i__ + j * b_dim1] = alpha * temp;
159 for (j = 1; j <= i__1; ++j) {
161 for (i__ = 1; i__ <= i__2; ++i__) {
162 temp = b[i__ + j * b_dim1];
164 temp *= a[i__ + i__ * a_dim1];
167 for (k = i__ + 1; k <= i__3; ++k) {
168 temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
170 b[i__ + j * b_dim1] = alpha * temp;
176 if (*transa=='N' || *transa=='n') {
179 for (j = n; j >= 1; --j) {
182 temp *= a[j + j * a_dim1];
185 for (i__ = 1; i__ <= i__1; ++i__) {
186 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
189 for (k = 1; k <= i__1; ++k) {
190 if (fabs(a[k + j * a_dim1])>GMX_DOUBLE_MIN) {
191 temp = alpha * a[k + j * a_dim1];
193 for (i__ = 1; i__ <= i__2; ++i__) {
194 b[i__ + j * b_dim1] += temp * b[i__ + k *
202 for (j = 1; j <= i__1; ++j) {
205 temp *= a[j + j * a_dim1];
208 for (i__ = 1; i__ <= i__2; ++i__) {
209 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
212 for (k = j + 1; k <= i__2; ++k) {
213 if (fabs(a[k + j * a_dim1])>GMX_DOUBLE_MIN) {
214 temp = alpha * a[k + j * a_dim1];
216 for (i__ = 1; i__ <= i__3; ++i__) {
217 b[i__ + j * b_dim1] += temp * b[i__ + k *
228 for (k = 1; k <= i__1; ++k) {
230 for (j = 1; j <= i__2; ++j) {
231 if (fabs(a[j + k * a_dim1])>GMX_DOUBLE_MIN) {
232 temp = alpha * a[j + k * a_dim1];
234 for (i__ = 1; i__ <= i__3; ++i__) {
235 b[i__ + j * b_dim1] += temp * b[i__ + k *
242 temp *= a[k + k * a_dim1];
244 if (fabs(temp-1.0)>GMX_DOUBLE_EPS) {
246 for (i__ = 1; i__ <= i__2; ++i__) {
247 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
252 for (k = n; k >= 1; --k) {
254 for (j = k + 1; j <= i__1; ++j) {
255 if (fabs(a[j + k * a_dim1])>GMX_DOUBLE_MIN) {
256 temp = alpha * a[j + k * a_dim1];
258 for (i__ = 1; i__ <= i__2; ++i__) {
259 b[i__ + j * b_dim1] += temp * b[i__ + k *
266 temp *= a[k + k * a_dim1];
268 if (fabs(temp-1.0)>GMX_DOUBLE_EPS) {
270 for (i__ = 1; i__ <= i__1; ++i__) {
271 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];