Merge release-5-0 into master
[alexxy/gromacs.git] / src / gromacs / linearalgebra / gmx_blas / strmm.c
1 #include <math.h>
2
3 #include "gromacs/utility/real.h"
4 #include "../gmx_blas.h"
5
6 void 
7 F77_FUNC(strmm,STRMM)(const char *side, 
8                       const char *uplo, 
9                       const char *transa, 
10                       const char *diag, 
11                       int *m__, 
12                       int *n__, 
13                       float *alpha__, 
14                       float *a, 
15                       int *lda__, 
16                       float *b, 
17                       int *ldb__)
18 {
19     int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
20     
21     int m = *m__;
22     int n = *n__;
23     int lda = *lda__;
24     int ldb = *ldb__;
25     float alpha = *alpha__;
26     
27     /* Local variables */
28     int i__, j, k, info;
29     float temp;
30     int lside;
31     int nrowa;
32     int upper;
33     int nounit;
34     a_dim1 = lda;
35     a_offset = 1 + a_dim1;
36     a -= a_offset;
37     b_dim1 = ldb;
38     b_offset = 1 + b_dim1;
39     b -= b_offset;
40
41     /* Function Body */
42     lside = (*side=='L' || *side=='l');
43     if (lside) {
44         nrowa = m;
45     } else {
46         nrowa = n;
47     }
48     nounit = (*diag=='N' || *diag=='n');
49     upper = (*uplo=='U' || *uplo=='u');
50
51     info = 0;
52
53     if (n == 0) {
54         return;
55     }
56     if (fabs(alpha)<GMX_FLOAT_MIN) {
57         i__1 = n;
58         for (j = 1; j <= i__1; ++j) {
59             i__2 = m;
60             for (i__ = 1; i__ <= i__2; ++i__) {
61                 b[i__ + j * b_dim1] = 0.;
62             }
63         }
64         return;
65     }
66     if (lside) {
67         if (*transa=='N' || *transa=='n') {
68             if (upper) {
69                 i__1 = n;
70                 for (j = 1; j <= i__1; ++j) {
71                     i__2 = m;
72                     for (k = 1; k <= i__2; ++k) {
73                         if ( fabs(b[k + j * b_dim1])>GMX_FLOAT_MIN) {
74                             temp = alpha * b[k + j * b_dim1];
75                             i__3 = k - 1;
76                             for (i__ = 1; i__ <= i__3; ++i__) {
77                                 b[i__ + j * b_dim1] += temp * a[i__ + k * 
78                                         a_dim1];
79                             }
80                             if (nounit) {
81                                 temp *= a[k + k * a_dim1];
82                             }
83                             b[k + j * b_dim1] = temp;
84                         }
85                     }
86                 }
87             } else {
88                 i__1 = n;
89                 for (j = 1; j <= i__1; ++j) {
90                     for (k = m; k >= 1; --k) {
91                         if (fabs(b[k + j * b_dim1])>GMX_FLOAT_MIN) {
92                             temp = alpha * b[k + j * b_dim1];
93                             b[k + j * b_dim1] = temp;
94                             if (nounit) {
95                                 b[k + j * b_dim1] *= a[k + k * a_dim1];
96                             }
97                             i__2 = m;
98                             for (i__ = k + 1; i__ <= i__2; ++i__) {
99                                 b[i__ + j * b_dim1] += temp * a[i__ + k * 
100                                         a_dim1];
101                             }
102                         }
103                     }
104                 }
105             }
106         } else {
107
108             if (upper) {
109                 i__1 = n;
110                 for (j = 1; j <= i__1; ++j) {
111                     for (i__ = m; i__ >= 1; --i__) {
112                         temp = b[i__ + j * b_dim1];
113                         if (nounit) {
114                             temp *= a[i__ + i__ * a_dim1];
115                         }
116                         i__2 = i__ - 1;
117                         for (k = 1; k <= i__2; ++k) {
118                             temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
119                         }
120                         b[i__ + j * b_dim1] = alpha * temp;
121                     }
122                 }
123             } else {
124                 i__1 = n;
125                 for (j = 1; j <= i__1; ++j) {
126                     i__2 = m;
127                     for (i__ = 1; i__ <= i__2; ++i__) {
128                         temp = b[i__ + j * b_dim1];
129                         if (nounit) {
130                             temp *= a[i__ + i__ * a_dim1];
131                         }
132                         i__3 = m;
133                         for (k = i__ + 1; k <= i__3; ++k) {
134                             temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
135                         }
136                         b[i__ + j * b_dim1] = alpha * temp;
137                     }
138                 }
139             }
140         }
141     } else {
142         if (*transa=='N' || *transa=='n') {
143
144             if (upper) {
145                 for (j = n; j >= 1; --j) {
146                     temp = alpha;
147                     if (nounit) {
148                         temp *= a[j + j * a_dim1];
149                     }
150                     i__1 = m;
151                     for (i__ = 1; i__ <= i__1; ++i__) {
152                         b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
153                     }
154                     i__1 = j - 1;
155                     for (k = 1; k <= i__1; ++k) {
156                         if ( fabs(a[k + j * a_dim1])>GMX_FLOAT_MIN) {
157                             temp = alpha * a[k + j * a_dim1];
158                             i__2 = m;
159                             for (i__ = 1; i__ <= i__2; ++i__) {
160                                 b[i__ + j * b_dim1] += temp * b[i__ + k * 
161                                         b_dim1];
162                             }
163                         }
164                     }
165                 }
166             } else {
167                 i__1 = n;
168                 for (j = 1; j <= i__1; ++j) {
169                     temp = alpha;
170                     if (nounit) {
171                         temp *= a[j + j * a_dim1];
172                     }
173                     i__2 = m;
174                     for (i__ = 1; i__ <= i__2; ++i__) {
175                         b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
176                     }
177                     i__2 = n;
178                     for (k = j + 1; k <= i__2; ++k) {
179                         if ( fabs(a[k + j * a_dim1])>GMX_FLOAT_MIN) {
180                             temp = alpha * a[k + j * a_dim1];
181                             i__3 = m;
182                             for (i__ = 1; i__ <= i__3; ++i__) {
183                                 b[i__ + j * b_dim1] += temp * b[i__ + k * 
184                                         b_dim1];
185                             }
186                         }
187                     }
188                 }
189             }
190         } else {
191
192             if (upper) {
193                 i__1 = n;
194                 for (k = 1; k <= i__1; ++k) {
195                     i__2 = k - 1;
196                     for (j = 1; j <= i__2; ++j) {
197                         if ( fabs(a[j + k * a_dim1])>GMX_FLOAT_MIN ) {
198                             temp = alpha * a[j + k * a_dim1];
199                             i__3 = m;
200                             for (i__ = 1; i__ <= i__3; ++i__) {
201                                 b[i__ + j * b_dim1] += temp * b[i__ + k * 
202                                         b_dim1];
203                             }
204                         }
205                     }
206                     temp = alpha;
207                     if (nounit) {
208                         temp *= a[k + k * a_dim1];
209                     }
210                     if ( fabs(temp-1.0)>GMX_FLOAT_EPS) {
211                         i__2 = m;
212                         for (i__ = 1; i__ <= i__2; ++i__) {
213                             b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
214                         }
215                     }
216                 }
217             } else {
218                 for (k = n; k >= 1; --k) {
219                     i__1 = n;
220                     for (j = k + 1; j <= i__1; ++j) {
221                         if ( fabs(a[j + k * a_dim1])>GMX_FLOAT_MIN) {
222                             temp = alpha * a[j + k * a_dim1];
223                             i__2 = m;
224                             for (i__ = 1; i__ <= i__2; ++i__) {
225                                 b[i__ + j * b_dim1] += temp * b[i__ + k * 
226                                         b_dim1];
227                             }
228                         }
229                     }
230                     temp = alpha;
231                     if (nounit) {
232                         temp *= a[k + k * a_dim1];
233                     }
234                     if ( fabs(temp-1.0)>GMX_FLOAT_EPS) {
235                         i__1 = m;
236                         for (i__ = 1; i__ <= i__1; ++i__) {
237                             b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
238                         }
239                     }
240                 }
241             }
242         }
243     }
244
245     return;
246
247 }
248
249