Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : These files are semi-automatic translations by f2c from the original netlib BLAS library.
3 : The source has been modified to (mostly) use modern C formatting, and to get rid of
4 : compiler warnings. Any errors in doing this should be blamed on the GROMACS developers, and
5 : not the reference BLAS implementation.
6 :
7 : The reference BLAS implementation is available from http://www.netlib.org/blas
8 :
9 : BLAS does not come with a formal named "license", but a general statement that
10 :
11 : "The reference BLAS is a freely-available software package. It is available from netlib
12 : via anonymous ftp and the World Wide Web. Thus, it can be included in commercial software
13 : packages (and has been). We only ask that proper credit be given to the authors."
14 :
15 : While the rest of GROMACS is LGPL, we think it's only fair to give you the same rights to
16 : our modified BLAS files as the original netlib versions, so do what you want with them.
17 : However, be warned that we have only tested that they to the right thing in the cases used
18 : in GROMACS (primarily full & sparse matrix diagonalization), so in most cases it is a much
19 : better idea to use the full reference implementation.
20 :
21 : Erik Lindahl, 2008-10-07.
22 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
23 : #if ! defined (__PLUMED_HAS_EXTERNAL_BLAS)
24 : #include <cmath>
25 : #include "blas.h"
26 :
27 : namespace PLMD{
28 : namespace blas{
29 : double
30 0 : PLUMED_BLAS_F77_FUNC(dasum,DASUM)(int *n__,
31 : double *dx,
32 : int *incx__)
33 : {
34 : int i__1, i__2;
35 :
36 : int i__, m, mp1;
37 : double dtemp;
38 : int nincx;
39 :
40 0 : int n = *n__;
41 0 : int incx = *incx__;
42 :
43 0 : --dx;
44 :
45 : dtemp = 0.;
46 0 : if (n <= 0 || incx <= 0) {
47 : return 0.0;
48 : }
49 0 : if (incx != 1) {
50 0 : nincx = n * incx;
51 : i__1 = nincx;
52 : i__2 = incx;
53 0 : for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
54 0 : dtemp += std::abs(dx[i__]);
55 : }
56 : return dtemp;
57 : }
58 :
59 0 : m = n % 6;
60 0 : if (m != 0) {
61 : i__2 = m;
62 0 : for (i__ = 1; i__ <= i__2; ++i__) {
63 0 : dtemp += std::abs(dx[i__]);
64 : }
65 0 : if (n < 6) {
66 : return dtemp;
67 : }
68 : }
69 0 : mp1 = m + 1;
70 : i__2 = n;
71 0 : for (i__ = mp1; i__ <= i__2; i__ += 6) {
72 0 : dtemp = dtemp + std::abs(dx[i__]) + std::abs(dx[i__ + 1]) +
73 0 : std::abs(dx[i__ + 2]) + std::abs(dx[i__+ 3]) + std::abs(dx[i__ + 4]) +
74 0 : std::abs(dx[i__ + 5]);
75 : }
76 : return dtemp;
77 : }
78 :
79 :
80 : }
81 : }
82 : #include "blas.h"
83 :
84 :
85 : namespace PLMD{
86 : namespace blas{
87 : void
88 1538741 : PLUMED_BLAS_F77_FUNC(daxpy,DAXPY)(int * n_arg,
89 : double * da_arg,
90 : double * dx,
91 : int * incx_arg,
92 : double * dy,
93 : int * incy_arg)
94 : {
95 : int i,ix,iy;
96 1538741 : int n=*n_arg;
97 1538741 : double da=*da_arg;
98 1538741 : int incx = *incx_arg;
99 1538741 : int incy = *incy_arg;
100 :
101 1538741 : if (n<=0)
102 : return;
103 :
104 1538741 : if(incx!=1 || incy!=1) {
105 : ix = 0;
106 : iy = 0;
107 0 : if(incx<0)
108 0 : ix = (1-n)*incx;
109 0 : if(incy<0)
110 0 : iy = (1-n)*incy;
111 :
112 0 : for(i=0;i<n;i++,ix+=incx,iy+=incy)
113 0 : dy[iy] += da*dx[ix];
114 :
115 : return;
116 :
117 : } else {
118 :
119 : /* unroll */
120 :
121 8639360 : for(i=0;i<(n-4);i+=4) {
122 7100619 : dy[i] += da*dx[i];
123 7100619 : dy[i+1] += da*dx[i+1];
124 7100619 : dy[i+2] += da*dx[i+2];
125 7100619 : dy[i+3] += da*dx[i+3];
126 : }
127 : /* continue with current value of i */
128 5393300 : for(;i<n;i++)
129 3854559 : dy[i] += da*dx[i];
130 : }
131 : }
132 : }
133 : }
134 : #include "blas.h"
135 :
136 : namespace PLMD{
137 : namespace blas{
138 : void
139 4471174 : PLUMED_BLAS_F77_FUNC(dcopy,DCOPY)(int *n__,
140 : double *dx,
141 : int *incx__,
142 : double *dy,
143 : int *incy__)
144 : {
145 : int i,ix,iy;
146 :
147 4471174 : int n= *n__;
148 4471174 : int incx = *incx__;
149 4471174 : int incy = *incy__;
150 :
151 :
152 4471174 : if(incx!=1 || incy!=1) {
153 : ix = 0;
154 : iy = 0;
155 5884 : if(incx<0)
156 0 : ix = (1-n)*(incx);
157 5884 : if(incy<0)
158 0 : iy = (1-n)*(incy);
159 :
160 1156404 : for(i=0;i<n;i++,ix+=incx,iy+=incy)
161 1150520 : dy[iy] = dx[ix];
162 :
163 : return;
164 :
165 : } else {
166 :
167 : /* unroll */
168 :
169 4572411 : for(i=0;i<(n-8);i+=8) {
170 107121 : dy[i] = dx[i];
171 107121 : dy[i+1] = dx[i+1];
172 107121 : dy[i+2] = dx[i+2];
173 107121 : dy[i+3] = dx[i+3];
174 107121 : dy[i+4] = dx[i+4];
175 107121 : dy[i+5] = dx[i+5];
176 107121 : dy[i+6] = dx[i+6];
177 107121 : dy[i+7] = dx[i+7];
178 : }
179 : /* continue with current value of i */
180 18533861 : for(;i<n;i++)
181 14068571 : dy[i] = dx[i];
182 : }
183 : }
184 : }
185 : }
186 : #include "blas.h"
187 :
188 : namespace PLMD{
189 : namespace blas{
190 : double
191 1249223 : PLUMED_BLAS_F77_FUNC(ddot,DDOT)(int *n_arg,
192 : double *dx,
193 : int *incx_arg,
194 : double *dy,
195 : int *incy_arg)
196 : {
197 : int i,ix,iy,m;
198 1249223 : int n=*n_arg;
199 1249223 : int incx = *incx_arg;
200 1249223 : int incy = *incy_arg;
201 : double t1;
202 :
203 1249223 : if(n<=0)
204 : return 0.0;
205 :
206 : t1 = 0.0;
207 :
208 1249223 : if(incx!=1 || incy!=1) {
209 : ix = 0;
210 : iy = 0;
211 0 : if(incx<0)
212 0 : ix = (1-n)*incx;
213 0 : if(incy<0)
214 0 : iy = (1-n)*incy;
215 :
216 0 : for(i=0;i<n;i++,ix+=incx,iy+=incy)
217 0 : t1 += dx[ix] * dy[iy];
218 :
219 : return t1;
220 :
221 : } else {
222 :
223 1249223 : m = n%5;
224 :
225 4371031 : for(i=0;i<m;i++)
226 3121808 : t1 += dx[i] * dy[i];
227 :
228 : /* unroll */
229 1285544 : for(i=m;i<n;i+=5)
230 36321 : t1 = t1 + dx[i] * dy[i]
231 36321 : + dx[i+1] * dy[i+1]
232 36321 : + dx[i+2] * dy[i+2]
233 36321 : + dx[i+3] * dy[i+3]
234 36321 : + dx[i+4] * dy[i+4];
235 :
236 : return t1;
237 : }
238 : }
239 :
240 :
241 : }
242 : }
243 : #include <cctype>
244 : #include <cmath>
245 :
246 : #include "real.h"
247 :
248 : #include "blas.h"
249 :
250 : namespace PLMD{
251 : namespace blas{
252 : void
253 504 : PLUMED_BLAS_F77_FUNC(dgemm,DGEMM)(const char *transa,
254 : const char *transb,
255 : int *m__,
256 : int *n__,
257 : int *k__,
258 : double *alpha__,
259 : double *a,
260 : int *lda__,
261 : double *b,
262 : int *ldb__,
263 : double *beta__,
264 : double *c,
265 : int *ldc__)
266 : {
267 504 : const char tra=std::toupper(*transa);
268 504 : const char trb=std::toupper(*transb);
269 : double temp;
270 : int i,j,l;
271 :
272 504 : int m = *m__;
273 504 : int n = *n__;
274 504 : int k = *k__;
275 504 : int lda = *lda__;
276 504 : int ldb = *ldb__;
277 504 : int ldc = *ldc__;
278 :
279 504 : double alpha = *alpha__;
280 504 : double beta = *beta__;
281 :
282 504 : if(m==0 || n==0 || (( std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN || k==0) && std::abs(beta-1.0)<PLUMED_GMX_DOUBLE_EPS))
283 : return;
284 :
285 455 : if(std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN) {
286 0 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN) {
287 0 : for(j=0;j<n;j++)
288 0 : for(i=0;i<m;i++)
289 0 : c[j*(ldc)+i] = 0.0;
290 : } else {
291 : /* nonzero beta */
292 0 : for(j=0;j<n;j++)
293 0 : for(i=0;i<m;i++)
294 0 : c[j*(ldc)+i] *= beta;
295 : }
296 0 : return;
297 : }
298 :
299 455 : if(trb=='N') {
300 346 : if(tra=='N') {
301 :
302 17396 : for(j=0;j<n;j++) {
303 17120 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN) {
304 934292 : for(i=0;i<m;i++)
305 924369 : c[j*(ldc)+i] = 0.0;
306 7197 : } else if(std::abs(beta-1.0)>PLUMED_GMX_DOUBLE_EPS) {
307 0 : for(i=0;i<m;i++)
308 0 : c[j*(ldc)+i] *= beta;
309 : }
310 898860 : for(l=0;l<k;l++) {
311 881740 : if( std::abs(b[ j*(ldb) + l ])>PLUMED_GMX_DOUBLE_MIN) {
312 879210 : temp = alpha * b[ j*(ldb) + l ];
313 223022936 : for(i=0;i<m;i++)
314 222143726 : c[j*(ldc)+i] += temp * a[l*(lda)+i];
315 : }
316 : }
317 : }
318 : } else {
319 : /* transpose A, but not B */
320 2098 : for(j=0;j<n;j++) {
321 579738 : for(i=0;i<m;i++) {
322 : temp = 0.0;
323 125807406 : for(l=0;l<k;l++)
324 125229696 : temp += a[i*(lda)+l] * b[j*(ldb)+l];
325 577710 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
326 0 : c[j*(ldc)+i] = alpha * temp;
327 : else
328 577710 : c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
329 : }
330 : }
331 : }
332 : } else {
333 : /* transpose B */
334 109 : if(tra=='N') {
335 :
336 : /* transpose B, but not A */
337 :
338 23285 : for(j=0;j<n;j++) {
339 23176 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN) {
340 0 : for(i=0;i<m;i++)
341 0 : c[j*(ldc)+i] = 0.0;
342 23176 : } else if(std::abs(beta-1.0)>PLUMED_GMX_DOUBLE_EPS) {
343 0 : for(i=0;i<m;i++)
344 0 : c[j*(ldc)+i] *= beta;
345 : }
346 831190 : for(l=0;l<k;l++) {
347 808014 : if( std::abs(b[ l*(ldb) + j ])>PLUMED_GMX_DOUBLE_MIN) {
348 701310 : temp = alpha * b[ l*(ldb) + j ];
349 207288274 : for(i=0;i<m;i++)
350 206586964 : c[j*(ldc)+i] += temp * a[l*(lda)+i];
351 : }
352 : }
353 : }
354 :
355 : } else {
356 : /* Transpose both A and B */
357 0 : for(j=0;j<n;j++) {
358 0 : for(i=0;i<m;i++) {
359 : temp = 0.0;
360 0 : for(l=0;l<k;l++)
361 0 : temp += a[i*(lda)+l] * b[l*(ldb)+j];
362 0 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
363 0 : c[j*(ldc)+i] = alpha * temp;
364 : else
365 0 : c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
366 : }
367 : }
368 : }
369 : }
370 : }
371 : }
372 : }
373 : #include <cctype>
374 : #include <cmath>
375 :
376 : #include "real.h"
377 :
378 : #include "blas.h"
379 :
380 : namespace PLMD{
381 : namespace blas{
382 : void
383 1287182 : PLUMED_BLAS_F77_FUNC(dgemv,DGEMV)(const char *trans,
384 : int *m__,
385 : int *n__,
386 : double *alpha__,
387 : double *a,
388 : int *lda__,
389 : double *x,
390 : int *incx__,
391 : double *beta__,
392 : double *y,
393 : int *incy__)
394 : {
395 1287182 : const char ch=std::toupper(*trans);
396 : int lenx,leny,kx,ky;
397 : int i,j,jx,jy,ix,iy;
398 : double temp;
399 :
400 1287182 : int m = *m__;
401 1287182 : int n = *n__;
402 1287182 : double alpha = *alpha__;
403 1287182 : double beta = *beta__;
404 1287182 : int incx = *incx__;
405 1287182 : int incy = *incy__;
406 1287182 : int lda = *lda__;
407 :
408 1287182 : if(n<=0 || m<=0 || (std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN && std::abs(beta-1.0)<PLUMED_GMX_DOUBLE_EPS))
409 : return;
410 :
411 1286990 : if(ch=='N') {
412 : lenx = n;
413 : leny = m;
414 : } else {
415 : lenx = m;
416 : leny = n;
417 : }
418 :
419 1286990 : if(incx>0)
420 : kx = 1;
421 : else
422 0 : kx = 1 - (lenx -1)*(incx);
423 :
424 1286990 : if(incy>0)
425 : ky = 1;
426 : else
427 0 : ky = 1 - (leny -1)*(incy);
428 :
429 1286990 : if(std::abs(beta-1.0)>PLUMED_GMX_DOUBLE_EPS) {
430 1265500 : if(incy==1) {
431 1265500 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
432 3333508 : for(i=0;i<leny;i++)
433 2068008 : y[i] = 0.0;
434 : else
435 0 : for(i=0;i<leny;i++)
436 0 : y[i] *= beta;
437 : } else {
438 : /* non-unit incr. */
439 : iy = ky;
440 0 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
441 0 : for(i=0;i<leny;i++,iy+=incy)
442 0 : y[iy] = 0.0;
443 : else
444 0 : for(i=0;i<leny;i++,iy+=incy)
445 0 : y[iy] *= beta;
446 : }
447 : }
448 :
449 1286990 : if(std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN)
450 : return;
451 :
452 1286990 : if(ch=='N') {
453 : jx = kx;
454 27222 : if(incy==1) {
455 511632 : for(j=1;j<=n;j++,jx+=incx)
456 484794 : if(std::abs(x[jx-1])>PLUMED_GMX_DOUBLE_MIN) {
457 484794 : temp = alpha * x[jx-1];
458 63100049 : for(i=1;i<=m;i++)
459 62615255 : y[i-1] += temp * a[(j-1)*(lda)+(i-1)];
460 : }
461 : } else {
462 : /* non-unit y incr. */
463 6720 : for(j=1;j<=n;j++,jx+=incx)
464 6336 : if(std::abs(x[jx-1])>PLUMED_GMX_DOUBLE_MIN) {
465 6336 : temp = alpha * x[jx-1];
466 : iy = ky;
467 1921920 : for(i=1;i<=m;i++,iy+=incy)
468 1915584 : y[iy-1] += temp * a[(j-1)*(lda)+(i-1)];
469 : }
470 : }
471 : } else {
472 : /* transpose */
473 : jy = ky;
474 1259768 : if(incx==1) {
475 3230172 : for(j=1;j<=n;j++,jy+=incy) {
476 : temp = 0.0;
477 62212273 : for(i=1;i<=m;i++)
478 60241113 : temp += a[(j-1)*(lda)+(i-1)] * x[i-1];
479 1971160 : y[jy-1] += alpha * temp;
480 : }
481 : } else {
482 : /* non-unit y incr. */
483 121296 : for(j=1;j<=n;j++,jy+=incy) {
484 : temp = 0.0;
485 : ix = kx;
486 3833628 : for(i=1;i<=m;i++,ix+=incx)
487 3713088 : temp += a[(j-1)*(lda)+(i-1)] * x[ix-1];
488 120540 : y[jy-1] += alpha * temp;
489 : }
490 : }
491 : }
492 : }
493 :
494 : }
495 : }
496 : #include <cmath>
497 :
498 : #include "real.h"
499 :
500 : #include "blas.h"
501 :
502 : namespace PLMD{
503 : namespace blas{
504 : void
505 1275614 : PLUMED_BLAS_F77_FUNC(dger,DGER)(int *m__,
506 : int *n__,
507 : double *alpha__,
508 : double *x,
509 : int *incx__,
510 : double *y,
511 : int *incy__,
512 : double *a,
513 : int *lda__)
514 : {
515 : int ix,kx,jy;
516 : int i,j;
517 : double temp;
518 :
519 :
520 1275614 : int m = *m__;
521 1275614 : int n = *n__;
522 1275614 : int incx = *incx__;
523 1275614 : int incy = *incy__;
524 1275614 : int lda = *lda__;
525 1275614 : double alpha = *alpha__;
526 :
527 1275614 : if(m<=0 || n<=0 || std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN)
528 : return;
529 :
530 1275614 : if(incy>0)
531 : jy = 0;
532 : else
533 0 : jy = incy * (1 - n);
534 :
535 1275614 : if(incx==1) {
536 3021858 : for(j=0;j<n;j++,jy+=incy)
537 1746244 : if(std::abs(y[jy])>PLUMED_GMX_DOUBLE_MIN) {
538 1746234 : temp = alpha * y[jy];
539 8618460 : for(i=0;i<m;i++)
540 6872226 : a[j*(lda)+i] += temp*x[i];
541 : }
542 : } else {
543 : /* non-unit incx */
544 0 : if(incx>0)
545 : kx = 0;
546 : else
547 0 : kx = incx * (1 - m);
548 :
549 0 : for(j=0;j<n;j++,jy+=incy) {
550 0 : if(std::abs(y[jy])>PLUMED_GMX_DOUBLE_MIN) {
551 0 : temp = alpha * y[jy];
552 : ix = kx;
553 0 : for(i=0;i<m;i++,ix+=incx)
554 0 : a[j*(lda)+i] += temp*x[ix];
555 : }
556 : }
557 : }
558 : return;
559 : }
560 : }
561 : }
562 : #include <cmath>
563 :
564 : #include "real.h"
565 : #include "blas.h"
566 :
567 : namespace PLMD{
568 : namespace blas{
569 : double
570 1262812 : PLUMED_BLAS_F77_FUNC(dnrm2,DNRM2)(int * n__,
571 : double * x,
572 : int * incx__)
573 : {
574 : int ix,max_ix;
575 : double ssq,scale,absxi,t;
576 :
577 1262812 : int n = *n__;
578 1262812 : int incx = *incx__;
579 :
580 1262812 : if(n<1 || incx<1)
581 : return 0;
582 1262812 : else if (n==1) {
583 625509 : t = x[0];
584 625509 : if(t>=0)
585 : return t;
586 : else
587 464697 : return -t;
588 : }
589 :
590 : scale = 0.0;
591 : ssq = 1.0;
592 :
593 637303 : max_ix = 1+(n-1)*(incx);
594 3211597 : for(ix=1;ix<=max_ix;ix+=incx) {
595 2574294 : t = x[ix-1];
596 2574294 : if(std::abs(t)>PLUMED_GMX_DOUBLE_MIN) {
597 2572105 : absxi = (t>=0) ? t : (-t);
598 2572105 : if(scale<absxi) {
599 922099 : t = scale/absxi;
600 922099 : t = t*t;
601 922099 : ssq = ssq*t + 1.0;
602 : scale = absxi;
603 : } else {
604 1650006 : t = absxi/scale;
605 1650006 : ssq += t*t;
606 : }
607 : }
608 : }
609 637303 : return scale*std::sqrt(ssq);
610 :
611 : }
612 :
613 :
614 :
615 : }
616 : }
617 : #include "blas.h"
618 :
619 : namespace PLMD{
620 : namespace blas{
621 : void
622 1800 : PLUMED_BLAS_F77_FUNC(drot,DROT)(int *n__,
623 : double *dx,
624 : int *incx__,
625 : double *dy,
626 : int *incy__,
627 : double *c__,
628 : double *s__)
629 : {
630 : int i,ix,iy;
631 : double dtemp;
632 :
633 1800 : int n = *n__;
634 1800 : int incx = *incx__;
635 1800 : int incy = *incy__;
636 1800 : double c = *c__;
637 1800 : double s = *s__;
638 :
639 1800 : if(incx!=1 || incy!=1) {
640 : ix = 0;
641 : iy = 0;
642 900 : if(incx<0)
643 0 : ix = (1-n)*(incx);
644 900 : if(incy<0)
645 0 : iy = (1-n)*(incy);
646 :
647 6016 : for(i=0;i<n;i++,ix+=incx,iy+=incy) {
648 5116 : dtemp = (c) * dx[ix] + (s) * dy[iy];
649 5116 : dy[iy] = (c) * dy[iy] - (s) * dx[ix];
650 5116 : dx[ix] = dtemp;
651 : }
652 :
653 : return;
654 :
655 : } else {
656 :
657 : /* unit increments */
658 5955 : for(i=0;i<n;i++) {
659 5055 : dtemp = (c) * dx[i] + (s) * dy[i];
660 5055 : dy[i] = (c) * dy[i] - (s) * dx[i];
661 5055 : dx[i] = dtemp;
662 : }
663 :
664 : }
665 : }
666 : }
667 : }
668 : #include "blas.h"
669 :
670 : namespace PLMD{
671 : namespace blas{
672 : void
673 2137688 : PLUMED_BLAS_F77_FUNC(dscal,DSCAL)(int * n__,
674 : double * fact__,
675 : double * dx,
676 : int * incx__)
677 : {
678 : int nincx,i;
679 :
680 2137688 : int n = *n__;
681 2137688 : double fact = *fact__;
682 2137688 : int incx = *incx__;
683 :
684 2137688 : if(n<=0 || incx<=0)
685 : return;
686 :
687 2121177 : if(incx==1) {
688 : /* Unrool factor 5 */
689 2303046 : for(i=0;i<(n-5);i+=5) {
690 185904 : dx[i] *= fact;
691 185904 : dx[i+1] *= fact;
692 185904 : dx[i+2] *= fact;
693 185904 : dx[i+3] *= fact;
694 185904 : dx[i+4] *= fact;
695 : }
696 : /* continue with current value of i */
697 7302624 : for(;i<n;i++)
698 5185482 : dx[i] *= fact;
699 :
700 : return;
701 : } else {
702 : /* inc != 1 */
703 4035 : nincx = n * (incx);
704 162438 : for (i=0;i<nincx;i+=incx)
705 158403 : dx[i] *= fact;
706 :
707 : return;
708 : }
709 :
710 : }
711 : }
712 : }
713 : #include "blas.h"
714 :
715 : namespace PLMD{
716 : namespace blas{
717 : void
718 8490 : PLUMED_BLAS_F77_FUNC(dswap,DSWAP)(int *n__,
719 : double *dx,
720 : int *incx__,
721 : double *dy,
722 : int *incy__)
723 : {
724 : int i,ix,iy;
725 : double d1,d2,d3;
726 :
727 8490 : int n = *n__;
728 8490 : int incx = *incx__;
729 8490 : int incy = *incy__;
730 :
731 8490 : if(n<=0)
732 : return;
733 :
734 8490 : if(incx==1 && incy==1) {
735 77201 : for(i=0;i<(n-3);i+=3) {
736 72777 : d1 = dx[i];
737 72777 : d2 = dx[i+1];
738 72777 : d3 = dx[i+2];
739 72777 : dx[i] = dy[i];
740 72777 : dx[i+1] = dy[i+1];
741 72777 : dx[i+2] = dy[i+2];
742 72777 : dy[i] = d1;
743 72777 : dy[i+1] = d2;
744 72777 : dy[i+2] = d3;
745 : }
746 : /* continue with last i value */
747 14349 : for(;i<n;i++) {
748 9925 : d1 = dx[i];
749 9925 : dx[i] = dy[i];
750 9925 : dy[i] = d1;
751 : }
752 :
753 : } else {
754 : ix = 0;
755 : iy = 0;
756 4066 : if(incx<0)
757 0 : ix = incx * (1 - n);
758 4066 : if(incy<0)
759 0 : iy = incy * (1 - n);
760 :
761 200881 : for(i=0;i<n;i++,ix+=incx,iy+=incy) {
762 196815 : d1 = dx[ix];
763 196815 : dx[ix] = dy[iy];
764 196815 : dy[iy] = d1;
765 : }
766 : }
767 : return;
768 : }
769 :
770 : }
771 : }
772 : #include <cctype>
773 : #include <cmath>
774 :
775 : #include "real.h"
776 : #include "blas.h"
777 :
778 : namespace PLMD{
779 : namespace blas{
780 : void
781 1249222 : PLUMED_BLAS_F77_FUNC(dsymv,DSYMV)(const char *uplo,
782 : int *n__,
783 : double *alpha__,
784 : double *a,
785 : int *lda__,
786 : double *x,
787 : int *incx__,
788 : double *beta__,
789 : double *y,
790 : int *incy__)
791 : {
792 1249222 : const char ch=std::toupper(*uplo);
793 : int kx,ky,i,j,ix,iy,jx,jy;
794 : double temp1,temp2;
795 :
796 1249222 : int n = *n__;
797 1249222 : int lda = *lda__;
798 1249222 : int incx = *incx__;
799 1249222 : int incy = *incy__;
800 1249222 : double alpha = *alpha__;
801 1249222 : double beta = *beta__;
802 :
803 1249222 : if(n<=0 || incx==0 || incy==0)
804 : return;
805 :
806 1249222 : if(incx>0)
807 : kx = 1;
808 : else
809 0 : kx = 1 - (n -1)*(incx);
810 :
811 1249222 : if(incy>0)
812 : ky = 1;
813 : else
814 0 : ky = 1 - (n -1)*(incy);
815 :
816 1249222 : if(std::abs(beta-1.0)>PLUMED_GMX_DOUBLE_EPS) {
817 1249222 : if(incy==1) {
818 1249222 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
819 4552630 : for(i=1;i<=n;i++)
820 3303408 : y[i-1] = 0.0;
821 : else
822 0 : for(i=1;i<=n;i++)
823 0 : y[i-1] *= beta;
824 : } else {
825 : /* non-unit incr. */
826 : iy = ky;
827 0 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
828 0 : for(i=1;i<=n;i++) {
829 0 : y[iy-1] = 0.0;
830 0 : iy += incy;
831 : }
832 : else
833 0 : for(i=1;i<=n;i++) {
834 0 : y[iy-1] *= beta;
835 0 : iy += incy;
836 : }
837 : }
838 : }
839 :
840 1249222 : if(std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN)
841 : return;
842 :
843 1249222 : if(ch=='U') {
844 1249222 : if(incx==1 && incy==1) {
845 4552630 : for(j=1;j<=n;j++) {
846 3303408 : temp1 = alpha * x[j-1];
847 : temp2 = 0.0;
848 29831600 : for(i=1;i<j;i++) {
849 26528192 : y[i-1] += temp1*a[(j-1)*(lda)+(i-1)];
850 26528192 : temp2 += a[(j-1)*(lda)+(i-1)] * x[i-1];
851 : }
852 3303408 : y[j-1] += temp1*a[(j-1)*(lda)+(j-1)] + alpha *temp2;
853 : }
854 : } else {
855 : /* non-unit incr. */
856 : jx = kx;
857 : jy = ky;
858 0 : for(j=1;j<=n;j++) {
859 0 : temp1 = alpha * x[jx-1];
860 : temp2 = 0.0;
861 : ix = kx;
862 : iy = ky;
863 0 : for(i=1;i<j;i++) {
864 0 : y[iy-1] += temp1 * a[(j-1)*(lda)+(i-1)];
865 0 : temp2 += a[(j-1)*(lda)+(i-1)] * x[ix-1];
866 0 : ix += incx;
867 0 : iy += incy;
868 : }
869 0 : y[jy-1] += temp1*a[(j-1)*(lda)+(j-1)] + alpha*temp2;
870 0 : jx += incx;
871 0 : jy += incy;
872 : }
873 : }
874 : } else {
875 : /* lower */
876 0 : if(incx==1 && incy==1) {
877 0 : for(j=1;j<=n;j++) {
878 0 : temp1 = alpha * x[j-1];
879 : temp2 = 0.0;
880 0 : y[j-1] += temp1 * a[(j-1)*(lda)+(j-1)];
881 0 : for(i=j+1;i<=n;i++) {
882 0 : y[i-1] += temp1*a[(j-1)*(lda)+(i-1)];
883 0 : temp2 += a[(j-1)*(lda)+(i-1)] * x[i-1];
884 : }
885 0 : y[j-1] += alpha *temp2;
886 : }
887 : } else {
888 : /* non-unit incr. */
889 : jx = kx;
890 : jy = ky;
891 0 : for(j=1;j<=n;j++) {
892 0 : temp1 = alpha * x[jx-1];
893 : temp2 = 0.0;
894 0 : y[jy-1] += temp1 * a[(j-1)*(lda)+(j-1)];
895 : ix = jx;
896 : iy = jy;
897 0 : for(i=j+1;i<=n;i++) {
898 0 : ix += incx;
899 0 : iy += incy;
900 0 : y[iy-1] += temp1 * a[(j-1)*(lda)+(i-1)];
901 0 : temp2 += a[(j-1)*(lda)+(i-1)] * x[ix-1];
902 : }
903 0 : y[jy-1] += alpha*temp2;
904 0 : jx += incx;
905 0 : jy += incy;
906 : }
907 : }
908 : }
909 : return;
910 : }
911 : }
912 : }
913 : #include <cctype>
914 : #include <cmath>
915 :
916 : #include "real.h"
917 :
918 : #include "blas.h"
919 :
920 : namespace PLMD{
921 : namespace blas{
922 : void
923 1248709 : PLUMED_BLAS_F77_FUNC(dsyr2,DSYR2)(const char * uplo,
924 : int * n__,
925 : double * alpha__,
926 : double * x,
927 : int * incx__,
928 : double * y,
929 : int * incy__,
930 : double * a,
931 : int * lda__)
932 : {
933 : int kx,ky,ix,iy,jx,jy,j,i;
934 : double temp1,temp2;
935 1248709 : const char ch=std::toupper(*uplo);
936 :
937 1248709 : int n = *n__;
938 1248709 : int lda = *lda__;
939 1248709 : int incx = *incx__;
940 1248709 : int incy = *incy__;
941 1248709 : float alpha = *alpha__;
942 :
943 :
944 1248709 : if(n<=0 || std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN || incx==0 || incy==0 ||
945 1248709 : (ch != 'U' && ch != 'L'))
946 : return;
947 :
948 : jx = jy = kx = ky = 0;
949 :
950 : /* init start points for non-unit increments */
951 1248709 : if(incx!=1 || incy!=1) {
952 0 : if(incx>0)
953 : kx = 1;
954 : else
955 0 : kx = 1 - (n - 1)*(incx);
956 0 : if(incy>0)
957 : ky = 1;
958 : else
959 0 : ky = 1 - (n - 1)*(incy);
960 :
961 : jx = kx;
962 : jy = ky;
963 : }
964 :
965 1248709 : if(ch == 'U') {
966 : /* Data in upper part of A */
967 1248709 : if(incx==1 && incy==1) {
968 : /* Unit increments for both x and y */
969 4410178 : for(j=1;j<=n;j++) {
970 3161469 : if( std::abs(x[j-1])>PLUMED_GMX_DOUBLE_MIN || std::abs(y[j-1])>PLUMED_GMX_DOUBLE_MIN ) {
971 3161467 : temp1 = alpha * y[j-1];
972 3161467 : temp2 = alpha * x[j-1];
973 10111771 : for(i=1;i<=j;i++)
974 6950304 : a[(j-1)*(lda)+(i-1)] += x[i-1]*temp1 + y[i-1]*temp2;
975 : }
976 : }
977 : } else {
978 :
979 : /* non-unit increments */
980 0 : for(j=1;j<=n;j++) {
981 :
982 0 : if( std::abs(x[jx-1])>PLUMED_GMX_DOUBLE_MIN || std::abs(y[jy-1])>PLUMED_GMX_DOUBLE_MIN ) {
983 0 : temp1 = alpha * y[jy-1];
984 0 : temp2 = alpha * x[jx-1];
985 : ix = kx;
986 : iy = ky;
987 0 : for(i=1;i<=j;i++) {
988 0 : a[(j-1)*(lda)+(i-1)] += x[ix-1]*temp1 + y[iy-1]*temp2;
989 0 : ix += incx;
990 0 : iy += incy;
991 : }
992 : }
993 0 : jx += incx;
994 0 : jy += incy;
995 : }
996 : }
997 : } else {
998 : /* Data in lower part of A */
999 0 : if(incx==1 && incy==1) {
1000 : /* Unit increments for both x and y */
1001 0 : for(j=1;j<=n;j++) {
1002 0 : if( std::abs(x[j-1])>PLUMED_GMX_DOUBLE_MIN || std::abs(y[j-1])>PLUMED_GMX_DOUBLE_MIN ) {
1003 0 : temp1 = alpha * y[j-1];
1004 0 : temp2 = alpha * x[j-1];
1005 0 : for(i=j;i<=n;i++)
1006 0 : a[(j-1)*(lda)+(i-1)] += x[i-1]*temp1 + y[i-1]*temp2;
1007 : }
1008 : }
1009 : } else {
1010 :
1011 : /* non-unit increments */
1012 0 : for(j=1;j<=n;j++) {
1013 :
1014 0 : if( std::abs(x[jx-1])>PLUMED_GMX_DOUBLE_MIN || std::abs(y[jy-1])>PLUMED_GMX_DOUBLE_MIN ) {
1015 0 : temp1 = alpha * y[jy-1];
1016 0 : temp2 = alpha * x[jx-1];
1017 : ix = jx;
1018 : iy = jy;
1019 0 : for(i=j;i<=n;i++) {
1020 0 : a[(j-1)*(lda)+(i-1)] += x[ix-1]*temp1 + y[iy-1]*temp2;
1021 0 : ix += incx;
1022 0 : iy += incy;
1023 : }
1024 : }
1025 0 : jx += incx;
1026 0 : jy += incy;
1027 : }
1028 : }
1029 : }
1030 :
1031 : return;
1032 : }
1033 : }
1034 : }
1035 : #include <cctype>
1036 : #include <cmath>
1037 :
1038 : #include "real.h"
1039 : #include "blas.h"
1040 :
1041 : namespace PLMD{
1042 : namespace blas{
1043 : void
1044 19 : PLUMED_BLAS_F77_FUNC(dsyr2k,DSYR2K)(const char *uplo,
1045 : const char *trans,
1046 : int *n__,
1047 : int *k__,
1048 : double *alpha__,
1049 : double *a,
1050 : int *lda__,
1051 : double *b,
1052 : int *ldb__,
1053 : double *beta__,
1054 : double *c,
1055 : int *ldc__)
1056 : {
1057 : char ch1,ch2;
1058 : int i,j,l;
1059 : double temp1,temp2;
1060 :
1061 :
1062 19 : int n = *n__;
1063 19 : int k = *k__;
1064 19 : int lda = *lda__;
1065 19 : int ldb = *ldb__;
1066 19 : int ldc = *ldc__;
1067 :
1068 19 : double alpha = *alpha__;
1069 19 : double beta = *beta__;
1070 :
1071 19 : ch1 = std::toupper(*uplo);
1072 19 : ch2 = std::toupper(*trans);
1073 :
1074 19 : if(n==0 || ( ( std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN || k==0 ) && std::abs(beta-1.0)<PLUMED_GMX_DOUBLE_EPS))
1075 : return;
1076 :
1077 19 : if(std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN ) {
1078 0 : if(ch1=='U') {
1079 0 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
1080 0 : for(j=1;j<=n;j++)
1081 0 : for(i=1;i<=j;i++)
1082 0 : c[(j-1)*(ldc)+(i-1)] = 0.0;
1083 : else
1084 0 : for(j=1;j<=n;j++)
1085 0 : for(i=1;i<=j;i++)
1086 0 : c[(j-1)*(ldc)+(i-1)] *= beta;
1087 : } else {
1088 : /* lower */
1089 0 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
1090 0 : for(j=1;j<=n;j++)
1091 0 : for(i=j;i<=n;i++)
1092 0 : c[(j-1)*(ldc)+(i-1)] = 0.0;
1093 : else
1094 0 : for(j=1;j<=n;j++)
1095 0 : for(i=j;i<=n;i++)
1096 0 : c[(j-1)*(ldc)+(i-1)] *= beta;
1097 : }
1098 0 : return;
1099 : }
1100 :
1101 19 : if(ch2=='N') {
1102 19 : if(ch1=='U') {
1103 5029 : for(j=1;j<=n;j++) {
1104 5010 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
1105 0 : for(i=1;i<=j;i++)
1106 0 : c[(j-1)*(ldc)+(i-1)] = 0.0;
1107 5010 : else if(std::abs(beta-1.0)>PLUMED_GMX_DOUBLE_EPS)
1108 0 : for(i=1;i<=j;i++)
1109 0 : c[(j-1)*(ldc)+(i-1)] *= beta;
1110 140280 : for(l=1;l<=k;l++) {
1111 135270 : if( std::abs(a[(l-1)*(lda)+(j-1)])>PLUMED_GMX_DOUBLE_MIN ||
1112 0 : std::abs(b[(l-1)*(ldb)+(j-1)])>PLUMED_GMX_DOUBLE_MIN) {
1113 135270 : temp1 = alpha * b[(l-1)*(ldb)+(j-1)];
1114 135270 : temp2 = alpha * a[(l-1)*(lda)+(j-1)];
1115 21195810 : for(i=1;i<=j;i++)
1116 21060540 : c[(j-1)*(ldc)+(i-1)] +=
1117 21060540 : a[(l-1)*(lda)+(i-1)] * temp1 +
1118 21060540 : b[(l-1)*(ldb)+(i-1)] * temp2;
1119 : }
1120 : }
1121 : }
1122 : } else {
1123 : /* lower */
1124 0 : for(j=1;j<=n;j++) {
1125 0 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
1126 0 : for(i=j;i<=n;i++)
1127 0 : c[(j-1)*(ldc)+(i-1)] = 0.0;
1128 0 : else if(std::abs(beta-1.0)>PLUMED_GMX_DOUBLE_EPS)
1129 0 : for(i=j;i<=n;i++)
1130 0 : c[(j-1)*(ldc)+(i-1)] *= beta;
1131 0 : for(l=1;l<=k;l++) {
1132 0 : if( std::abs(a[(l-1)*(lda)+(j-1)])>PLUMED_GMX_DOUBLE_MIN ||
1133 0 : std::abs(b[(l-1)*(ldb)+(j-1)])>PLUMED_GMX_DOUBLE_MIN) {
1134 0 : temp1 = alpha * b[(l-1)*(ldb)+(j-1)];
1135 0 : temp2 = alpha * a[(l-1)*(lda)+(j-1)];
1136 0 : for(i=j;i<=n;i++)
1137 0 : c[(j-1)*(ldc)+(i-1)] +=
1138 0 : a[(l-1)*(lda)+(i-1)] * temp1 +
1139 0 : b[(l-1)*(ldb)+(i-1)] * temp2;
1140 : }
1141 : }
1142 : }
1143 : }
1144 : } else {
1145 : /* transpose */
1146 0 : if(ch1=='U') {
1147 0 : for(j=1;j<=n;j++)
1148 0 : for(i=1;i<=j;i++) {
1149 : temp1 = 0.0;
1150 : temp2 = 0.0;
1151 0 : for (l=1;l<=k;l++) {
1152 0 : temp1 += a[(i-1)*(lda)+(l-1)] * b[(j-1)*(ldb)+(l-1)];
1153 0 : temp2 += b[(i-1)*(ldb)+(l-1)] * a[(j-1)*(lda)+(l-1)];
1154 : }
1155 0 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
1156 0 : c[(j-1)*(ldc)+(i-1)] = alpha * (temp1 + temp2);
1157 : else
1158 0 : c[(j-1)*(ldc)+(i-1)] = beta * c[(j-1)*(ldc)+(i-1)] +
1159 0 : alpha * (temp1 + temp2);
1160 : }
1161 : } else {
1162 : /* lower */
1163 0 : for(j=1;j<=n;j++)
1164 0 : for(i=j;i<=n;i++) {
1165 : temp1 = 0.0;
1166 : temp2 = 0.0;
1167 0 : for (l=1;l<=k;l++) {
1168 0 : temp1 += a[(i-1)*(lda)+(l-1)] * b[(j-1)*(ldb)+(l-1)];
1169 0 : temp2 += b[(i-1)*(ldb)+(l-1)] * a[(j-1)*(lda)+(l-1)];
1170 : }
1171 0 : if(std::abs(beta)<PLUMED_GMX_DOUBLE_MIN)
1172 0 : c[(j-1)*(ldc)+(i-1)] = alpha * (temp1 + temp2);
1173 : else
1174 0 : c[(j-1)*(ldc)+(i-1)] = beta * c[(j-1)*(ldc)+(i-1)] +
1175 0 : alpha * (temp1 + temp2);
1176 : }
1177 : }
1178 : }
1179 : return;
1180 : }
1181 : }
1182 : }
1183 : #include <cmath>
1184 :
1185 : #include "real.h"
1186 :
1187 : #include "blas.h"
1188 :
1189 : namespace PLMD{
1190 : namespace blas{
1191 : void
1192 411 : PLUMED_BLAS_F77_FUNC(dtrmm,DTRMM)(const char *side,
1193 : const char *uplo,
1194 : const char *transa,
1195 : const char *diag,
1196 : int *m__,
1197 : int *n__,
1198 : double *alpha__,
1199 : double *a,
1200 : int *lda__,
1201 : double *b,
1202 : int *ldb__)
1203 : {
1204 : int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
1205 :
1206 411 : int m = *m__;
1207 411 : int n = *n__;
1208 411 : int lda = *lda__;
1209 411 : int ldb = *ldb__;
1210 411 : double alpha = *alpha__;
1211 :
1212 : /* Local variables */
1213 : int i__, j, k;
1214 : double temp;
1215 : int lside;
1216 : int upper;
1217 : int nounit;
1218 : a_dim1 = lda;
1219 411 : a_offset = 1 + a_dim1;
1220 411 : a -= a_offset;
1221 : b_dim1 = ldb;
1222 411 : b_offset = 1 + b_dim1;
1223 411 : b -= b_offset;
1224 :
1225 : /* Function Body */
1226 411 : lside = (*side=='L' || *side=='l');
1227 :
1228 411 : nounit = (*diag=='N' || *diag=='n');
1229 411 : upper = (*uplo=='U' || *uplo=='u');
1230 :
1231 411 : if (n == 0) {
1232 : return;
1233 : }
1234 411 : if (std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN) {
1235 : i__1 = n;
1236 0 : for (j = 1; j <= i__1; ++j) {
1237 : i__2 = m;
1238 0 : for (i__ = 1; i__ <= i__2; ++i__) {
1239 0 : b[i__ + j * b_dim1] = 0.;
1240 : }
1241 : }
1242 : return;
1243 : }
1244 411 : if (lside) {
1245 0 : if (*transa=='N' || *transa=='n') {
1246 0 : if (upper) {
1247 : i__1 = n;
1248 0 : for (j = 1; j <= i__1; ++j) {
1249 : i__2 = m;
1250 0 : for (k = 1; k <= i__2; ++k) {
1251 0 : if (std::abs(b[k + j * b_dim1])>PLUMED_GMX_DOUBLE_MIN) {
1252 0 : temp = alpha * b[k + j * b_dim1];
1253 : i__3 = k - 1;
1254 0 : for (i__ = 1; i__ <= i__3; ++i__) {
1255 0 : b[i__ + j * b_dim1] += temp * a[i__ + k * a_dim1];
1256 : }
1257 0 : if (nounit) {
1258 0 : temp *= a[k + k * a_dim1];
1259 : }
1260 0 : b[k + j * b_dim1] = temp;
1261 : }
1262 : }
1263 : }
1264 : } else {
1265 : i__1 = n;
1266 0 : for (j = 1; j <= i__1; ++j) {
1267 0 : for (k = m; k >= 1; --k) {
1268 0 : if (std::abs(b[k + j * b_dim1])>PLUMED_GMX_DOUBLE_MIN) {
1269 0 : temp = alpha * b[k + j * b_dim1];
1270 0 : b[k + j * b_dim1] = temp;
1271 0 : if (nounit) {
1272 0 : b[k + j * b_dim1] *= a[k + k * a_dim1];
1273 : }
1274 : i__2 = m;
1275 0 : for (i__ = k + 1; i__ <= i__2; ++i__) {
1276 0 : b[i__ + j * b_dim1] += temp * a[i__ + k *
1277 0 : a_dim1];
1278 : }
1279 : }
1280 : }
1281 : }
1282 : }
1283 : } else {
1284 :
1285 0 : if (upper) {
1286 : i__1 = n;
1287 0 : for (j = 1; j <= i__1; ++j) {
1288 0 : for (i__ = m; i__ >= 1; --i__) {
1289 0 : temp = b[i__ + j * b_dim1];
1290 0 : if (nounit) {
1291 0 : temp *= a[i__ + i__ * a_dim1];
1292 : }
1293 : i__2 = i__ - 1;
1294 0 : for (k = 1; k <= i__2; ++k) {
1295 0 : temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
1296 : }
1297 0 : b[i__ + j * b_dim1] = alpha * temp;
1298 : }
1299 : }
1300 : } else {
1301 : i__1 = n;
1302 0 : for (j = 1; j <= i__1; ++j) {
1303 : i__2 = m;
1304 0 : for (i__ = 1; i__ <= i__2; ++i__) {
1305 0 : temp = b[i__ + j * b_dim1];
1306 0 : if (nounit) {
1307 0 : temp *= a[i__ + i__ * a_dim1];
1308 : }
1309 : i__3 = m;
1310 0 : for (k = i__ + 1; k <= i__3; ++k) {
1311 0 : temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
1312 : }
1313 0 : b[i__ + j * b_dim1] = alpha * temp;
1314 : }
1315 : }
1316 : }
1317 : }
1318 : } else {
1319 411 : if (*transa=='N' || *transa=='n') {
1320 :
1321 137 : if (upper) {
1322 2376 : for (j = n; j >= 1; --j) {
1323 : temp = alpha;
1324 2287 : if (nounit) {
1325 0 : temp *= a[j + j * a_dim1];
1326 : }
1327 : i__1 = m;
1328 630793 : for (i__ = 1; i__ <= i__1; ++i__) {
1329 628506 : b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
1330 : }
1331 : i__1 = j - 1;
1332 36173 : for (k = 1; k <= i__1; ++k) {
1333 33886 : if (std::abs(a[k + j * a_dim1])>PLUMED_GMX_DOUBLE_MIN) {
1334 33254 : temp = alpha * a[k + j * a_dim1];
1335 : i__2 = m;
1336 9535918 : for (i__ = 1; i__ <= i__2; ++i__) {
1337 9502664 : b[i__ + j * b_dim1] += temp * b[i__ + k *
1338 9502664 : b_dim1];
1339 : }
1340 : }
1341 : }
1342 : }
1343 : } else {
1344 : i__1 = n;
1345 1124 : for (j = 1; j <= i__1; ++j) {
1346 : temp = alpha;
1347 1076 : if (nounit) {
1348 0 : temp *= a[j + j * a_dim1];
1349 : }
1350 : i__2 = m;
1351 271892 : for (i__ = 1; i__ <= i__2; ++i__) {
1352 270816 : b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
1353 : }
1354 : i__2 = n;
1355 16778 : for (k = j + 1; k <= i__2; ++k) {
1356 15702 : if (std::abs(a[k + j * a_dim1])>PLUMED_GMX_DOUBLE_MIN) {
1357 15486 : temp = alpha * a[k + j * a_dim1];
1358 : i__3 = m;
1359 4113622 : for (i__ = 1; i__ <= i__3; ++i__) {
1360 4098136 : b[i__ + j * b_dim1] += temp * b[i__ + k *
1361 4098136 : b_dim1];
1362 : }
1363 : }
1364 : }
1365 : }
1366 : }
1367 : } else {
1368 :
1369 274 : if (upper) {
1370 : i__1 = n;
1371 4471 : for (k = 1; k <= i__1; ++k) {
1372 : i__2 = k - 1;
1373 67533 : for (j = 1; j <= i__2; ++j) {
1374 63239 : if (std::abs(a[j + k * a_dim1])>PLUMED_GMX_DOUBLE_MIN) {
1375 62486 : temp = alpha * a[j + k * a_dim1];
1376 : i__3 = m;
1377 17678442 : for (i__ = 1; i__ <= i__3; ++i__) {
1378 17615956 : b[i__ + j * b_dim1] += temp * b[i__ + k *
1379 17615956 : b_dim1];
1380 : }
1381 : }
1382 : }
1383 : temp = alpha;
1384 4294 : if (nounit) {
1385 2007 : temp *= a[k + k * a_dim1];
1386 : }
1387 4294 : if (std::abs(temp-1.0)>PLUMED_GMX_DOUBLE_EPS) {
1388 : i__2 = m;
1389 538339 : for (i__ = 1; i__ <= i__2; ++i__) {
1390 536332 : b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
1391 : }
1392 : }
1393 : }
1394 : } else {
1395 2529 : for (k = n; k >= 1; --k) {
1396 : i__1 = n;
1397 38369 : for (j = k + 1; j <= i__1; ++j) {
1398 35937 : if (std::abs(a[j + k * a_dim1])>PLUMED_GMX_DOUBLE_MIN) {
1399 34810 : temp = alpha * a[j + k * a_dim1];
1400 : i__2 = m;
1401 9565722 : for (i__ = 1; i__ <= i__2; ++i__) {
1402 9530912 : b[i__ + j * b_dim1] += temp * b[i__ + k *
1403 9530912 : b_dim1];
1404 : }
1405 : }
1406 : }
1407 : temp = alpha;
1408 2432 : if (nounit) {
1409 1356 : temp *= a[k + k * a_dim1];
1410 : }
1411 2432 : if (std::abs(temp-1.0)>PLUMED_GMX_DOUBLE_EPS) {
1412 : i__1 = m;
1413 364346 : for (i__ = 1; i__ <= i__1; ++i__) {
1414 362990 : b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
1415 : }
1416 : }
1417 : }
1418 : }
1419 : }
1420 : }
1421 :
1422 : return;
1423 :
1424 : }
1425 :
1426 :
1427 : }
1428 : }
1429 : #include <cmath>
1430 :
1431 : #include "real.h"
1432 : #include "blas.h"
1433 :
1434 : namespace PLMD{
1435 : namespace blas{
1436 : void
1437 36253 : PLUMED_BLAS_F77_FUNC(dtrmv,DTRMV)(const char *uplo,
1438 : const char *trans,
1439 : const char *diag,
1440 : int *n__,
1441 : double *a,
1442 : int *lda__,
1443 : double *x,
1444 : int *incx__)
1445 : {
1446 : int a_dim1, a_offset, i__1, i__2;
1447 :
1448 : int i__, j, ix, jx, kx;
1449 : double temp;
1450 : int nounit;
1451 :
1452 36253 : int n = *n__;
1453 36253 : int lda = *lda__;
1454 36253 : int incx = *incx__;
1455 :
1456 : a_dim1 = lda;
1457 36253 : a_offset = 1 + a_dim1;
1458 36253 : a -= a_offset;
1459 36253 : --x;
1460 :
1461 36253 : if (n == 0) {
1462 : return;
1463 : }
1464 :
1465 19658 : nounit = (*diag=='n' || *diag=='N');
1466 :
1467 19658 : if (incx <= 0) {
1468 0 : kx = 1 - (n - 1) * incx;
1469 : } else {
1470 : kx = 1;
1471 : }
1472 :
1473 19658 : if (*trans=='N' || *trans=='n') {
1474 :
1475 19658 : if (*uplo=='U' || *uplo=='u') {
1476 18407 : if (incx == 1) {
1477 : i__1 = n;
1478 64159 : for (j = 1; j <= i__1; ++j) {
1479 45752 : if (std::abs(x[j])>PLUMED_GMX_DOUBLE_MIN) {
1480 : temp = x[j];
1481 : i__2 = j - 1;
1482 335151 : for (i__ = 1; i__ <= i__2; ++i__) {
1483 289615 : x[i__] += temp * a[i__ + j * a_dim1];
1484 : }
1485 45536 : if (nounit) {
1486 45536 : x[j] *= a[j + j * a_dim1];
1487 : }
1488 : }
1489 : }
1490 : } else {
1491 : jx = kx;
1492 : i__1 = n;
1493 0 : for (j = 1; j <= i__1; ++j) {
1494 0 : if (std::abs(x[jx])>PLUMED_GMX_DOUBLE_MIN) {
1495 : temp = x[jx];
1496 : ix = kx;
1497 : i__2 = j - 1;
1498 0 : for (i__ = 1; i__ <= i__2; ++i__) {
1499 0 : x[ix] += temp * a[i__ + j * a_dim1];
1500 0 : ix += incx;
1501 : }
1502 0 : if (nounit) {
1503 0 : x[jx] *= a[j + j * a_dim1];
1504 : }
1505 : }
1506 0 : jx += incx;
1507 : }
1508 : }
1509 : } else {
1510 1251 : if (incx == 1) {
1511 20575 : for (j = n; j >= 1; --j) {
1512 19324 : if (std::abs(x[j])>PLUMED_GMX_DOUBLE_MIN) {
1513 : temp = x[j];
1514 : i__1 = j + 1;
1515 208252 : for (i__ = n; i__ >= i__1; --i__) {
1516 188928 : x[i__] += temp * a[i__ + j * a_dim1];
1517 : }
1518 19324 : if (nounit) {
1519 19324 : x[j] *= a[j + j * a_dim1];
1520 : }
1521 : }
1522 : }
1523 : } else {
1524 0 : kx += (n - 1) * incx;
1525 : jx = kx;
1526 0 : for (j = n; j >= 1; --j) {
1527 0 : if (std::abs(x[jx])>PLUMED_GMX_DOUBLE_MIN) {
1528 : temp = x[jx];
1529 : ix = kx;
1530 : i__1 = j + 1;
1531 0 : for (i__ = n; i__ >= i__1; --i__) {
1532 0 : x[ix] += temp * a[i__ + j * a_dim1];
1533 0 : ix -= incx;
1534 : }
1535 0 : if (nounit) {
1536 0 : x[jx] *= a[j + j * a_dim1];
1537 : }
1538 : }
1539 0 : jx -= incx;
1540 : }
1541 : }
1542 : }
1543 : } else {
1544 :
1545 0 : if (*uplo=='U' || *uplo=='u') {
1546 0 : if (incx == 1) {
1547 0 : for (j = n; j >= 1; --j) {
1548 0 : temp = x[j];
1549 0 : if (nounit) {
1550 0 : temp *= a[j + j * a_dim1];
1551 : }
1552 0 : for (i__ = j - 1; i__ >= 1; --i__) {
1553 0 : temp += a[i__ + j * a_dim1] * x[i__];
1554 : }
1555 0 : x[j] = temp;
1556 : }
1557 : } else {
1558 0 : jx = kx + (n - 1) * incx;
1559 0 : for (j = n; j >= 1; --j) {
1560 0 : temp = x[jx];
1561 : ix = jx;
1562 0 : if (nounit) {
1563 0 : temp *= a[j + j * a_dim1];
1564 : }
1565 0 : for (i__ = j - 1; i__ >= 1; --i__) {
1566 0 : ix -= incx;
1567 0 : temp += a[i__ + j * a_dim1] * x[ix];
1568 : }
1569 0 : x[jx] = temp;
1570 0 : jx -= incx;
1571 : }
1572 : }
1573 : } else {
1574 0 : if (incx == 1) {
1575 : i__1 = n;
1576 0 : for (j = 1; j <= i__1; ++j) {
1577 0 : temp = x[j];
1578 0 : if (nounit) {
1579 0 : temp *= a[j + j * a_dim1];
1580 : }
1581 : i__2 = n;
1582 0 : for (i__ = j + 1; i__ <= i__2; ++i__) {
1583 0 : temp += a[i__ + j * a_dim1] * x[i__];
1584 : }
1585 0 : x[j] = temp;
1586 : }
1587 : } else {
1588 : jx = kx;
1589 : i__1 = n;
1590 0 : for (j = 1; j <= i__1; ++j) {
1591 0 : temp = x[jx];
1592 : ix = jx;
1593 0 : if (nounit) {
1594 0 : temp *= a[j + j * a_dim1];
1595 : }
1596 : i__2 = n;
1597 0 : for (i__ = j + 1; i__ <= i__2; ++i__) {
1598 0 : ix += incx;
1599 0 : temp += a[i__ + j * a_dim1] * x[ix];
1600 : }
1601 0 : x[jx] = temp;
1602 0 : jx += incx;
1603 : }
1604 : }
1605 : }
1606 : }
1607 :
1608 : return;
1609 :
1610 : }
1611 :
1612 :
1613 : }
1614 : }
1615 : #include <cctype>
1616 : #include <cmath>
1617 :
1618 : #include "real.h"
1619 : #include "blas.h"
1620 :
1621 : namespace PLMD{
1622 : namespace blas{
1623 : void
1624 0 : PLUMED_BLAS_F77_FUNC(dtrsm,DTRSM)(const char * side,
1625 : const char * uplo,
1626 : const char * transa,
1627 : const char * diag,
1628 : int * m__,
1629 : int * n__,
1630 : double *alpha__,
1631 : double *a,
1632 : int * lda__,
1633 : double *b,
1634 : int * ldb__)
1635 : {
1636 0 : const char xside = std::toupper(*side);
1637 0 : const char xuplo = std::toupper(*uplo);
1638 0 : const char xtrans = std::toupper(*transa);
1639 0 : const char xdiag = std::toupper(*diag);
1640 : int i,j,k;
1641 : double temp;
1642 :
1643 :
1644 0 : int m = *m__;
1645 0 : int n = *n__;
1646 0 : int lda = *lda__;
1647 0 : int ldb = *ldb__;
1648 0 : double alpha = *alpha__;
1649 :
1650 0 : if(n<=0)
1651 : return;
1652 :
1653 0 : if(std::abs(alpha)<PLUMED_GMX_DOUBLE_MIN) {
1654 0 : for(j=0;j<n;j++)
1655 0 : for(i=0;i<m;i++)
1656 0 : b[j*(ldb)+i] = 0.0;
1657 : return;
1658 : }
1659 :
1660 0 : if(xside=='L') {
1661 : /* left side */
1662 0 : if(xtrans=='N') {
1663 : /* No transpose */
1664 0 : if(xuplo=='U') {
1665 : /* upper */
1666 0 : for(j=0;j<n;j++) {
1667 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_DOUBLE_EPS) {
1668 0 : for(i=0;i<m;i++)
1669 0 : b[j*(ldb)+i] *= alpha;
1670 : }
1671 0 : for(k=m-1;k>=0;k--) {
1672 0 : if(std::abs(b[j*(ldb)+k])>PLUMED_GMX_DOUBLE_MIN) {
1673 0 : if(xdiag=='N')
1674 0 : b[j*(ldb)+k] /= a[k*(lda)+k];
1675 0 : for(i=0;i<k;i++)
1676 0 : b[j*(ldb)+i] -= b[j*(ldb)+k]*a[k*(lda)+i];
1677 : }
1678 : }
1679 : }
1680 : } else {
1681 : /* lower */
1682 0 : for(j=0;j<n;j++) {
1683 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_DOUBLE_EPS)
1684 0 : for(i=0;i<m;i++)
1685 0 : b[j*(ldb)+i] *= alpha;
1686 0 : for(k=0;k<m;k++) {
1687 0 : if(std::abs(b[j*(ldb)+k])>PLUMED_GMX_DOUBLE_MIN) {
1688 0 : if(xdiag=='N')
1689 0 : b[j*(ldb)+k] /= a[k*(lda)+k];
1690 0 : for(i=k+1;i<m;i++)
1691 0 : b[j*(ldb)+i] -= b[j*(ldb)+k]*a[k*(lda)+i];
1692 : }
1693 : }
1694 : }
1695 : }
1696 : } else {
1697 : /* Transpose */
1698 0 : if(xuplo=='U') {
1699 : /* upper */
1700 0 : for(j=0;j<n;j++) {
1701 0 : for(i=0;i<m;i++) {
1702 0 : temp = alpha * b[j*(ldb)+i];
1703 0 : for(k=0;k<i;k++)
1704 0 : temp -= a[i*(lda)+k] * b[j*(ldb)+k];
1705 0 : if(xdiag=='N')
1706 0 : temp /= a[i*(lda)+i];
1707 0 : b[j*(ldb)+i] = temp;
1708 : }
1709 : }
1710 : } else {
1711 : /* lower */
1712 0 : for(j=0;j<n;j++) {
1713 0 : for(i=m-1;i>=0;i--) {
1714 0 : temp = alpha * b[j*(ldb)+i];
1715 0 : for(k=i+1;k<m;k++)
1716 0 : temp -= a[i*(lda)+k] * b[j*(ldb)+k];
1717 0 : if(xdiag=='N')
1718 0 : temp /= a[i*(lda)+i];
1719 0 : b[j*(ldb)+i] = temp;
1720 : }
1721 : }
1722 : }
1723 : }
1724 : } else {
1725 : /* right side */
1726 0 : if(xtrans=='N') {
1727 : /* No transpose */
1728 0 : if(xuplo=='U') {
1729 : /* upper */
1730 0 : for(j=0;j<n;j++) {
1731 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_DOUBLE_EPS)
1732 0 : for(i=0;i<m;i++)
1733 0 : b[j*(ldb)+i] *= alpha;
1734 0 : for(k=0;k<j;k++) {
1735 0 : if(std::abs(a[j*(lda)+k])>PLUMED_GMX_DOUBLE_MIN) {
1736 0 : for(i=0;i<m;i++)
1737 0 : b[j*(ldb)+i] -= a[j*(lda)+k]*b[k*(ldb)+i];
1738 : }
1739 : }
1740 0 : if(xdiag=='N') {
1741 0 : temp = 1.0/a[j*(lda)+j];
1742 0 : for(i=0;i<m;i++)
1743 0 : b[j*(ldb)+i] *= temp;
1744 : }
1745 : }
1746 : } else {
1747 : /* lower */
1748 0 : for(j=n-1;j>=0;j--) {
1749 0 : if(std::abs(alpha)>PLUMED_GMX_DOUBLE_MIN)
1750 0 : for(i=0;i<m;i++)
1751 0 : b[j*(ldb)+i] *= alpha;
1752 0 : for(k=j+1;k<n;k++) {
1753 0 : if(std::abs(a[j*(lda)+k])>PLUMED_GMX_DOUBLE_MIN) {
1754 0 : for(i=0;i<m;i++)
1755 0 : b[j*(ldb)+i] -= a[j*(lda)+k]*b[k*(ldb)+i];
1756 : }
1757 : }
1758 0 : if(xdiag=='N') {
1759 0 : temp = 1.0/a[j*(lda)+j];
1760 0 : for(i=0;i<m;i++)
1761 0 : b[j*(ldb)+i] *= temp;
1762 : }
1763 : }
1764 : }
1765 : } else {
1766 : /* Transpose */
1767 0 : if(xuplo=='U') {
1768 : /* upper */
1769 0 : for(k=n-1;k>=0;k--) {
1770 0 : if(xdiag=='N') {
1771 0 : temp = 1.0/a[k*(lda)+k];
1772 0 : for(i=0;i<m;i++)
1773 0 : b[k*(ldb)+i] *= temp;
1774 : }
1775 0 : for(j=0;j<k;j++) {
1776 0 : if(std::abs(a[k*(lda)+j])>PLUMED_GMX_DOUBLE_MIN) {
1777 : temp = a[k*(lda)+j];
1778 0 : for(i=0;i<m;i++)
1779 0 : b[j*(ldb)+i] -= temp * b[k*(ldb)+i];
1780 : }
1781 : }
1782 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_DOUBLE_EPS)
1783 0 : for(i=0;i<m;i++)
1784 0 : b[k*(ldb)+i] *= alpha;
1785 : }
1786 : } else {
1787 : /* lower */
1788 0 : for(k=0;k<n;k++) {
1789 0 : if(xdiag=='N') {
1790 0 : temp = 1.0/a[k*(lda)+k];
1791 0 : for(i=0;i<m;i++)
1792 0 : b[k*(ldb)+i] *= temp;
1793 : }
1794 0 : for(j=k+1;j<n;j++) {
1795 0 : if(std::abs(a[k*(lda)+j])>PLUMED_GMX_DOUBLE_MIN) {
1796 : temp = a[k*(lda)+j];
1797 0 : for(i=0;i<m;i++)
1798 0 : b[j*(ldb)+i] -= temp * b[k*(ldb)+i];
1799 : }
1800 : }
1801 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_DOUBLE_EPS)
1802 0 : for(i=0;i<m;i++)
1803 0 : b[k*(ldb)+i] *= alpha;
1804 : }
1805 : }
1806 : }
1807 : }
1808 : }
1809 : }
1810 : }
1811 : #include <cmath>
1812 : #include "blas.h"
1813 :
1814 : namespace PLMD{
1815 : namespace blas{
1816 : int
1817 33025 : PLUMED_BLAS_F77_FUNC(idamax,IDAMAX)(int *n__,
1818 : double *dx,
1819 : int *incx__)
1820 : {
1821 : int i,ix,idxmax;
1822 : double dmax,tmp;
1823 :
1824 33025 : int n = *n__;
1825 33025 : int incx = *incx__;
1826 :
1827 33025 : if(n<1 || incx<=0)
1828 : return -1;
1829 :
1830 33025 : if(n==1)
1831 : return 1;
1832 :
1833 16514 : dmax = std::abs(dx[0]);
1834 : idxmax = 1;
1835 :
1836 16514 : if(incx==1) {
1837 33034 : for(i=1;i<n;i++) {
1838 16520 : tmp = std::abs(dx[i]);
1839 16520 : if(tmp>dmax) {
1840 : dmax = tmp;
1841 6 : idxmax = i+1;
1842 : }
1843 : }
1844 : } else {
1845 : /* Non-unit increments */
1846 : ix = incx; /* this is really 0 + an increment */
1847 0 : for(i=1;i<n;i++,ix+=incx) {
1848 0 : tmp = std::abs(dx[ix]);
1849 0 : if(tmp>dmax) {
1850 : dmax = tmp;
1851 0 : idxmax = ix+1;
1852 : }
1853 : }
1854 : }
1855 : return idxmax;
1856 : }
1857 : }
1858 : }
1859 : #include <cmath>
1860 : #include "blas.h"
1861 :
1862 : namespace PLMD{
1863 : namespace blas{
1864 : int
1865 0 : PLUMED_BLAS_F77_FUNC(isamax,ISAMAX)(int *n__,
1866 : float *dx,
1867 : int *incx__)
1868 : {
1869 : int i,ix,idxmax;
1870 : float dmax,tmp;
1871 :
1872 0 : int n = *n__;
1873 0 : int incx = *incx__;
1874 :
1875 0 : if(n<1 || incx<=0)
1876 : return -1;
1877 :
1878 0 : if(n==1)
1879 : return 1;
1880 :
1881 0 : dmax = std::abs(dx[0]);
1882 : idxmax = 1;
1883 :
1884 0 : if(incx==1) {
1885 0 : for(i=1;i<n;i++) {
1886 0 : tmp = std::abs(dx[i]);
1887 0 : if(tmp>dmax) {
1888 : dmax = tmp;
1889 0 : idxmax = i+1;
1890 : }
1891 : }
1892 : } else {
1893 : /* Non-unit increments */
1894 : ix = incx; /* this is really 0 + an increment */
1895 0 : for(i=1;i<n;i++,ix+=incx) {
1896 0 : tmp = std::abs(dx[ix]);
1897 0 : if(tmp>dmax) {
1898 : dmax = tmp;
1899 0 : idxmax = ix+1;
1900 : }
1901 : }
1902 : }
1903 : return idxmax;
1904 : }
1905 : }
1906 : }
1907 : #include <cmath>
1908 : #include "blas.h"
1909 :
1910 : namespace PLMD{
1911 : namespace blas{
1912 : float
1913 0 : PLUMED_BLAS_F77_FUNC(sasum,SASUM)(int *n__,
1914 : float *dx,
1915 : int *incx__)
1916 : {
1917 : int i__1, i__2;
1918 :
1919 : int i__, m, mp1;
1920 : float dtemp;
1921 : int nincx;
1922 :
1923 0 : int n = *n__;
1924 0 : int incx = *incx__;
1925 :
1926 :
1927 0 : --dx;
1928 :
1929 : dtemp = 0.;
1930 0 : if (n <= 0 || incx <= 0) {
1931 : return 0.0;
1932 : }
1933 0 : if (incx != 1) {
1934 0 : nincx = n * incx;
1935 : i__1 = nincx;
1936 : i__2 = incx;
1937 0 : for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
1938 0 : dtemp += std::abs(dx[i__]);
1939 : }
1940 : return dtemp;
1941 : }
1942 :
1943 0 : m = n % 6;
1944 0 : if (m != 0) {
1945 : i__2 = m;
1946 0 : for (i__ = 1; i__ <= i__2; ++i__) {
1947 0 : dtemp += std::abs(dx[i__]);
1948 : }
1949 0 : if (n < 6) {
1950 : return dtemp;
1951 : }
1952 : }
1953 0 : mp1 = m + 1;
1954 : i__2 = n;
1955 0 : for (i__ = mp1; i__ <= i__2; i__ += 6) {
1956 0 : dtemp = dtemp + std::abs(dx[i__]) + std::abs(dx[i__ + 1]) +
1957 0 : std::abs(dx[i__ + 2]) + std::abs(dx[i__+ 3]) + std::abs(dx[i__ + 4]) +
1958 0 : std::abs(dx[i__ + 5]);
1959 : }
1960 : return dtemp;
1961 : }
1962 :
1963 :
1964 : }
1965 : }
1966 : #include "blas.h"
1967 :
1968 :
1969 : namespace PLMD{
1970 : namespace blas{
1971 : void
1972 0 : PLUMED_BLAS_F77_FUNC(saxpy,SAXPY)(int * n_arg,
1973 : float * da_arg,
1974 : float * dx,
1975 : int * incx_arg,
1976 : float * dy,
1977 : int * incy_arg)
1978 : {
1979 : int i,ix,iy;
1980 0 : int n=*n_arg;
1981 0 : float da=*da_arg;
1982 0 : int incx = *incx_arg;
1983 0 : int incy = *incy_arg;
1984 :
1985 0 : if (n<=0)
1986 : return;
1987 :
1988 0 : if(incx!=1 || incy!=1) {
1989 : ix = 0;
1990 : iy = 0;
1991 0 : if(incx<0)
1992 0 : ix = (1-n)*incx;
1993 0 : if(incy<0)
1994 0 : iy = (1-n)*incy;
1995 :
1996 0 : for(i=0;i<n;i++,ix+=incx,iy+=incy)
1997 0 : dy[iy] += da*dx[ix];
1998 :
1999 : return;
2000 :
2001 : } else {
2002 :
2003 : /* unroll */
2004 :
2005 0 : for(i=0;i<(n-4);i+=4) {
2006 0 : dy[i] += da*dx[i];
2007 0 : dy[i+1] += da*dx[i+1];
2008 0 : dy[i+2] += da*dx[i+2];
2009 0 : dy[i+3] += da*dx[i+3];
2010 : }
2011 : /* continue with current value of i */
2012 0 : for(;i<n;i++)
2013 0 : dy[i] += da*dx[i];
2014 : }
2015 : }
2016 : }
2017 : }
2018 : #include "blas.h"
2019 :
2020 : namespace PLMD{
2021 : namespace blas{
2022 : void
2023 0 : PLUMED_BLAS_F77_FUNC(scopy,SCOPY)(int *n__,
2024 : float *dx,
2025 : int *incx__,
2026 : float *dy,
2027 : int *incy__)
2028 : {
2029 : int i,ix,iy;
2030 :
2031 0 : int n= *n__;
2032 0 : int incx = *incx__;
2033 0 : int incy = *incy__;
2034 :
2035 0 : if(incx!=1 || incy!=1)
2036 : {
2037 : ix = 0;
2038 : iy = 0;
2039 0 : if(incx<0)
2040 0 : ix = (1-n)*(incx);
2041 0 : if(incy<0)
2042 0 : iy = (1-n)*(incy);
2043 :
2044 0 : for(i=0;i<n;i++,ix+=incx,iy+=incy)
2045 0 : dy[iy] = dx[ix];
2046 :
2047 : return;
2048 :
2049 : } else {
2050 :
2051 : /* unroll */
2052 :
2053 0 : for(i=0;i<(n-8);i+=8) {
2054 0 : dy[i] = dx[i];
2055 0 : dy[i+1] = dx[i+1];
2056 0 : dy[i+2] = dx[i+2];
2057 0 : dy[i+3] = dx[i+3];
2058 0 : dy[i+4] = dx[i+4];
2059 0 : dy[i+5] = dx[i+5];
2060 0 : dy[i+6] = dx[i+6];
2061 0 : dy[i+7] = dx[i+7];
2062 : }
2063 : /* continue with current value of i */
2064 0 : for(;i<n;i++)
2065 0 : dy[i] = dx[i];
2066 : }
2067 : }
2068 : }
2069 : }
2070 : #include "blas.h"
2071 :
2072 :
2073 : namespace PLMD{
2074 : namespace blas{
2075 : float
2076 1 : PLUMED_BLAS_F77_FUNC(sdot,SDOT)(int *n_arg,
2077 : float *dx,
2078 : int *incx_arg,
2079 : float *dy,
2080 : int *incy_arg)
2081 : {
2082 : int i,ix,iy,m;
2083 1 : int n=*n_arg;
2084 1 : int incx = *incx_arg;
2085 1 : int incy = *incy_arg;
2086 : float t1;
2087 :
2088 1 : if(n<=0)
2089 : return 0.0;
2090 :
2091 : t1 = 0.0;
2092 :
2093 1 : if(incx!=1 || incy!=1) {
2094 : ix = 0;
2095 : iy = 0;
2096 0 : if(incx<0)
2097 0 : ix = (1-n)*incx;
2098 0 : if(incy<0)
2099 0 : iy = (1-n)*incy;
2100 :
2101 0 : for(i=0;i<n;i++,ix+=incx,iy+=incy)
2102 0 : t1 += dx[ix] * dy[iy];
2103 :
2104 : return t1;
2105 :
2106 : } else {
2107 :
2108 1 : m = n%5;
2109 :
2110 1 : for(i=0;i<m;i++)
2111 0 : t1 += dx[i] * dy[i];
2112 :
2113 : /* unroll */
2114 2 : for(i=m;i<n;i+=5)
2115 1 : t1 = t1 + dx[i] * dy[i]
2116 1 : + dx[i+1] * dy[i+1]
2117 1 : + dx[i+2] * dy[i+2]
2118 1 : + dx[i+3] * dy[i+3]
2119 1 : + dx[i+4] * dy[i+4];
2120 :
2121 : return t1;
2122 : }
2123 : }
2124 :
2125 :
2126 : }
2127 : }
2128 : #include <cctype>
2129 : #include <cmath>
2130 :
2131 : #include "real.h"
2132 : #include "blas.h"
2133 :
2134 : namespace PLMD{
2135 : namespace blas{
2136 : void
2137 0 : PLUMED_BLAS_F77_FUNC(sgemm,SGEMM)(const char *transa,
2138 : const char *transb,
2139 : int *m__,
2140 : int *n__,
2141 : int *k__,
2142 : float *alpha__,
2143 : float *a,
2144 : int *lda__,
2145 : float *b,
2146 : int *ldb__,
2147 : float *beta__,
2148 : float *c,
2149 : int *ldc__)
2150 : {
2151 0 : const char tra=std::toupper(*transa);
2152 0 : const char trb=std::toupper(*transb);
2153 : float temp;
2154 : int i,j,l;
2155 :
2156 0 : int m = *m__;
2157 0 : int n = *n__;
2158 0 : int k = *k__;
2159 0 : int lda = *lda__;
2160 0 : int ldb = *ldb__;
2161 0 : int ldc = *ldc__;
2162 :
2163 0 : float alpha = *alpha__;
2164 0 : float beta = *beta__;
2165 :
2166 0 : if(m==0 || n==0 || (( std::abs(alpha)<PLUMED_GMX_FLOAT_MIN || k==0) && std::abs(beta-1.0)<PLUMED_GMX_FLOAT_EPS))
2167 : return;
2168 :
2169 0 : if(std::abs(alpha)<PLUMED_GMX_FLOAT_MIN) {
2170 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN) {
2171 0 : for(j=0;j<n;j++)
2172 0 : for(i=0;i<m;i++)
2173 0 : c[j*(ldc)+i] = 0.0;
2174 : } else {
2175 : /* nonzero beta */
2176 0 : for(j=0;j<n;j++)
2177 0 : for(i=0;i<m;i++)
2178 0 : c[j*(ldc)+i] *= beta;
2179 : }
2180 0 : return;
2181 : }
2182 :
2183 0 : if(trb=='N') {
2184 0 : if(tra=='N') {
2185 :
2186 0 : for(j=0;j<n;j++) {
2187 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN) {
2188 0 : for(i=0;i<m;i++)
2189 0 : c[j*(ldc)+i] = 0.0;
2190 0 : } else if(std::abs(beta-1.0)>PLUMED_GMX_FLOAT_EPS) {
2191 0 : for(i=0;i<m;i++)
2192 0 : c[j*(ldc)+i] *= beta;
2193 : }
2194 0 : for(l=0;l<k;l++) {
2195 0 : if( std::abs(b[ j*(ldb) + l ])>PLUMED_GMX_FLOAT_MIN) {
2196 0 : temp = alpha * b[ j*(ldb) + l ];
2197 0 : for(i=0;i<m;i++)
2198 0 : c[j*(ldc)+i] += temp * a[l*(lda)+i];
2199 : }
2200 : }
2201 : }
2202 : } else {
2203 : /* transpose A, but not B */
2204 0 : for(j=0;j<n;j++) {
2205 0 : for(i=0;i<m;i++) {
2206 : temp = 0.0;
2207 0 : for(l=0;l<k;l++)
2208 0 : temp += a[i*(lda)+l] * b[j*(ldb)+l];
2209 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
2210 0 : c[j*(ldc)+i] = alpha * temp;
2211 : else
2212 0 : c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
2213 : }
2214 : }
2215 : }
2216 : } else {
2217 : /* transpose B */
2218 0 : if(tra=='N') {
2219 :
2220 : /* transpose B, but not A */
2221 :
2222 0 : for(j=0;j<n;j++) {
2223 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN) {
2224 0 : for(i=0;i<m;i++)
2225 0 : c[j*(ldc)+i] = 0.0;
2226 0 : } else if(std::abs(beta-1.0)>PLUMED_GMX_FLOAT_EPS) {
2227 0 : for(i=0;i<m;i++)
2228 0 : c[j*(ldc)+i] *= beta;
2229 : }
2230 0 : for(l=0;l<k;l++) {
2231 0 : if( std::abs(b[ l*(ldb) + j ])>PLUMED_GMX_FLOAT_MIN) {
2232 0 : temp = alpha * b[ l*(ldb) + j ];
2233 0 : for(i=0;i<m;i++)
2234 0 : c[j*(ldc)+i] += temp * a[l*(lda)+i];
2235 : }
2236 : }
2237 : }
2238 :
2239 : } else {
2240 : /* Transpose both A and B */
2241 0 : for(j=0;j<n;j++) {
2242 0 : for(i=0;i<m;i++) {
2243 : temp = 0.0;
2244 0 : for(l=0;l<k;l++)
2245 0 : temp += a[i*(lda)+l] * b[l*(ldb)+j];
2246 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
2247 0 : c[j*(ldc)+i] = alpha * temp;
2248 : else
2249 0 : c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
2250 : }
2251 : }
2252 : }
2253 : }
2254 : }
2255 : }
2256 : }
2257 : #include <cctype>
2258 : #include <cmath>
2259 :
2260 : #include "real.h"
2261 : #include "blas.h"
2262 :
2263 : namespace PLMD{
2264 : namespace blas{
2265 : void
2266 0 : PLUMED_BLAS_F77_FUNC(sgemv,SGEMV)(const char *trans,
2267 : int *m__,
2268 : int *n__,
2269 : float *alpha__,
2270 : float *a,
2271 : int *lda__,
2272 : float *x,
2273 : int *incx__,
2274 : float *beta__,
2275 : float *y,
2276 : int *incy__)
2277 : {
2278 0 : const char ch=std::toupper(*trans);
2279 : int lenx,leny,kx,ky;
2280 : int i,j,jx,jy,ix,iy;
2281 : float temp;
2282 :
2283 0 : int m = *m__;
2284 0 : int n = *n__;
2285 0 : float alpha = *alpha__;
2286 0 : float beta = *beta__;
2287 0 : int incx = *incx__;
2288 0 : int incy = *incy__;
2289 0 : int lda = *lda__;
2290 :
2291 0 : if(n<=0 || m<=0 || (std::abs(alpha)<PLUMED_GMX_FLOAT_MIN && std::abs(beta-1.0)<PLUMED_GMX_FLOAT_EPS))
2292 : return;
2293 :
2294 0 : if(ch=='N') {
2295 : lenx = n;
2296 : leny = m;
2297 : } else {
2298 : lenx = m;
2299 : leny = n;
2300 : }
2301 :
2302 0 : if(incx>0)
2303 : kx = 1;
2304 : else
2305 0 : kx = 1 - (lenx -1)*(incx);
2306 :
2307 0 : if(incy>0)
2308 : ky = 1;
2309 : else
2310 0 : ky = 1 - (leny -1)*(incy);
2311 :
2312 0 : if(std::abs(beta-1.0)>PLUMED_GMX_FLOAT_EPS) {
2313 0 : if(incy==1) {
2314 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
2315 0 : for(i=0;i<leny;i++)
2316 0 : y[i] = 0.0;
2317 : else
2318 0 : for(i=0;i<leny;i++)
2319 0 : y[i] *= beta;
2320 : } else {
2321 : /* non-unit incr. */
2322 : iy = ky;
2323 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
2324 0 : for(i=0;i<leny;i++,iy+=incy)
2325 0 : y[iy] = 0.0;
2326 : else
2327 0 : for(i=0;i<leny;i++,iy+=incy)
2328 0 : y[iy] *= beta;
2329 : }
2330 : }
2331 :
2332 0 : if(std::abs(alpha)<PLUMED_GMX_FLOAT_MIN)
2333 : return;
2334 :
2335 0 : if(ch=='N') {
2336 : jx = kx;
2337 0 : if(incy==1) {
2338 0 : for(j=1;j<=n;j++,jx+=incx)
2339 0 : if( std::abs(x[jx-1])>PLUMED_GMX_FLOAT_MIN) {
2340 0 : temp = alpha * x[jx-1];
2341 0 : for(i=1;i<=m;i++)
2342 0 : y[i-1] += temp * a[(j-1)*(lda)+(i-1)];
2343 : }
2344 : } else {
2345 : /* non-unit y incr. */
2346 0 : for(j=1;j<=n;j++,jx+=incx)
2347 0 : if( std::abs(x[jx-1])>PLUMED_GMX_FLOAT_MIN) {
2348 0 : temp = alpha * x[jx-1];
2349 : iy = ky;
2350 0 : for(i=1;i<=m;i++,iy+=incy)
2351 0 : y[iy-1] += temp * a[(j-1)*(lda)+(i-1)];
2352 : }
2353 : }
2354 : } else {
2355 : /* transpose */
2356 : jy = ky;
2357 0 : if(incx==1) {
2358 0 : for(j=1;j<=n;j++,jy+=incy) {
2359 : temp = 0.0;
2360 0 : for(i=1;i<=m;i++)
2361 0 : temp += a[(j-1)*(lda)+(i-1)] * x[i-1];
2362 0 : y[jy-1] += alpha * temp;
2363 : }
2364 : } else {
2365 : /* non-unit y incr. */
2366 0 : for(j=1;j<=n;j++,jy+=incy) {
2367 : temp = 0.0;
2368 : ix = kx;
2369 0 : for(i=1;i<=m;i++,ix+=incx)
2370 0 : temp += a[(j-1)*(lda)+(i-1)] * x[ix-1];
2371 0 : y[jy-1] += alpha * temp;
2372 : }
2373 : }
2374 : }
2375 : }
2376 :
2377 : }
2378 : }
2379 : #include <cmath>
2380 :
2381 : #include "real.h"
2382 : #include "blas.h"
2383 :
2384 : namespace PLMD{
2385 : namespace blas{
2386 : void
2387 0 : PLUMED_BLAS_F77_FUNC(sger,SGER)(int *m__,
2388 : int *n__,
2389 : float *alpha__,
2390 : float *x,
2391 : int *incx__,
2392 : float *y,
2393 : int *incy__,
2394 : float *a,
2395 : int *lda__)
2396 : {
2397 : int ix,kx,jy;
2398 : int i,j;
2399 : float temp;
2400 :
2401 0 : int m = *m__;
2402 0 : int n = *n__;
2403 0 : int incx = *incx__;
2404 0 : int incy = *incy__;
2405 0 : int lda = *lda__;
2406 0 : float alpha = *alpha__;
2407 :
2408 0 : if(m<=0 || n<=0 || std::abs(alpha)<PLUMED_GMX_FLOAT_MIN)
2409 : return;
2410 :
2411 0 : if(incy>0)
2412 : jy = 0;
2413 : else
2414 0 : jy = incy * (1 - n);
2415 :
2416 0 : if(incx==1) {
2417 0 : for(j=0;j<n;j++,jy+=incy)
2418 0 : if(std::abs(y[jy])>PLUMED_GMX_FLOAT_MIN) {
2419 0 : temp = alpha * y[jy];
2420 0 : for(i=0;i<m;i++)
2421 0 : a[j*(lda)+i] += temp*x[i];
2422 : }
2423 : } else {
2424 : /* non-unit incx */
2425 0 : if(incx>0)
2426 : kx = 0;
2427 : else
2428 0 : kx = incx * (1 - m);
2429 :
2430 0 : for(j=0;j<n;j++,jy+=incy) {
2431 0 : if(std::abs(y[jy])>PLUMED_GMX_FLOAT_MIN) {
2432 0 : temp = alpha * y[jy];
2433 : ix = kx;
2434 0 : for(i=0;i<m;i++,ix+=incx)
2435 0 : a[j*(lda)+i] += temp*x[ix];
2436 : }
2437 : }
2438 : }
2439 : return;
2440 : }
2441 : }
2442 : }
2443 : #include <cmath>
2444 :
2445 :
2446 : #include "real.h"
2447 : #include "blas.h"
2448 :
2449 : namespace PLMD{
2450 : namespace blas{
2451 : float
2452 0 : PLUMED_BLAS_F77_FUNC(snrm2,SNRM2)(int * n__,
2453 : float * x,
2454 : int * incx__)
2455 : {
2456 : int ix,max_ix;
2457 : float ssq,scale,absxi,t;
2458 :
2459 0 : int n = *n__;
2460 0 : int incx = *incx__;
2461 :
2462 0 : if(n<1 || incx<1)
2463 : return 0;
2464 0 : else if (n==1) {
2465 0 : t = x[0];
2466 0 : if(t>=0)
2467 : return t;
2468 : else
2469 0 : return -t;
2470 : }
2471 :
2472 : scale = 0.0;
2473 : ssq = 1.0;
2474 :
2475 0 : max_ix = 1+(n-1)*(incx);
2476 0 : for(ix=1;ix<=max_ix;ix+=incx) {
2477 0 : t = x[ix-1];
2478 0 : if(std::abs(t)>PLUMED_GMX_FLOAT_MIN) {
2479 0 : absxi = (t>=0) ? t : (-t);
2480 0 : if(scale<absxi) {
2481 0 : t = scale/absxi;
2482 0 : t = t*t;
2483 0 : ssq = ssq*t + 1.0;
2484 : scale = absxi;
2485 : } else {
2486 0 : t = absxi/scale;
2487 0 : ssq += t*t;
2488 : }
2489 : }
2490 : }
2491 0 : return scale*std::sqrt(ssq);
2492 :
2493 : }
2494 :
2495 :
2496 :
2497 : }
2498 : }
2499 : #include "blas.h"
2500 :
2501 : namespace PLMD{
2502 : namespace blas{
2503 : void
2504 0 : PLUMED_BLAS_F77_FUNC(srot,SROT)(int *n__,
2505 : float *dx,
2506 : int *incx__,
2507 : float *dy,
2508 : int *incy__,
2509 : float *c__,
2510 : float *s__)
2511 : {
2512 : int i,ix,iy;
2513 : float dtemp;
2514 :
2515 0 : int n = *n__;
2516 0 : int incx = *incx__;
2517 0 : int incy = *incy__;
2518 0 : float c = *c__;
2519 0 : float s = *s__;
2520 :
2521 0 : if(incx!=1 || incy!=1) {
2522 : ix = 0;
2523 : iy = 0;
2524 0 : if(incx<0)
2525 0 : ix = (1-n)*(incx);
2526 0 : if(incy<0)
2527 0 : iy = (1-n)*(incy);
2528 :
2529 0 : for(i=0;i<n;i++,ix+=incx,iy+=incy) {
2530 0 : dtemp = (c) * dx[ix] + (s) * dy[iy];
2531 0 : dy[iy] = (c) * dy[iy] - (s) * dx[ix];
2532 0 : dx[ix] = dtemp;
2533 : }
2534 :
2535 : return;
2536 :
2537 : } else {
2538 :
2539 : /* unit increments */
2540 0 : for(i=0;i<n;i++) {
2541 0 : dtemp = (c) * dx[i] + (s) * dy[i];
2542 0 : dy[i] = (c) * dy[i] - (s) * dx[i];
2543 0 : dx[i] = dtemp;
2544 : }
2545 :
2546 : }
2547 : }
2548 : }
2549 : }
2550 : #include "blas.h"
2551 :
2552 : namespace PLMD{
2553 : namespace blas{
2554 : void
2555 0 : PLUMED_BLAS_F77_FUNC(sscal,SSCAL)(int *n__,
2556 : float *fact__,
2557 : float *dx,
2558 : int *incx__)
2559 : {
2560 : int nincx,i;
2561 :
2562 0 : int n = *n__;
2563 0 : float fact = *fact__;
2564 0 : int incx = *incx__;
2565 :
2566 0 : if(n<=0 || incx<=0)
2567 : return;
2568 :
2569 0 : if(incx==1) {
2570 : /* Unrool factor 5 */
2571 0 : for(i=0;i<(n-5);i+=5) {
2572 0 : dx[i] *= fact;
2573 0 : dx[i+1] *= fact;
2574 0 : dx[i+2] *= fact;
2575 0 : dx[i+3] *= fact;
2576 0 : dx[i+4] *= fact;
2577 : }
2578 : /* continue with current value of i */
2579 0 : for(;i<n;i++)
2580 0 : dx[i] *= fact;
2581 :
2582 : return;
2583 : } else {
2584 : /* inc != 1 */
2585 0 : nincx = n * (incx);
2586 0 : for (i=0;i<nincx;i+=incx)
2587 0 : dx[i] *= fact;
2588 :
2589 : return;
2590 : }
2591 :
2592 : }
2593 : }
2594 : }
2595 : #include "blas.h"
2596 :
2597 : namespace PLMD{
2598 : namespace blas{
2599 : void
2600 0 : PLUMED_BLAS_F77_FUNC(sswap,SSWAP)(int *n__,
2601 : float *dx,
2602 : int *incx__,
2603 : float *dy,
2604 : int *incy__)
2605 : {
2606 : int i,ix,iy;
2607 : float d1,d2,d3;
2608 :
2609 0 : int n = *n__;
2610 0 : int incx = *incx__;
2611 0 : int incy = *incy__;
2612 :
2613 0 : if(n<=0)
2614 : return;
2615 :
2616 0 : if(incx==1 && incy==1) {
2617 0 : for(i=0;i<(n-3);i+=3) {
2618 0 : d1 = dx[i];
2619 0 : d2 = dx[i+1];
2620 0 : d3 = dx[i+2];
2621 0 : dx[i] = dy[i];
2622 0 : dx[i+1] = dy[i+1];
2623 0 : dx[i+2] = dy[i+2];
2624 0 : dy[i] = d1;
2625 0 : dy[i+1] = d2;
2626 0 : dy[i+2] = d3;
2627 : }
2628 : /* continue with last i value */
2629 0 : for(;i<n;i++) {
2630 0 : d1 = dx[i];
2631 0 : dx[i] = dy[i];
2632 0 : dy[i] = d1;
2633 : }
2634 :
2635 : } else {
2636 : ix = 0;
2637 : iy = 0;
2638 0 : if(incx<0)
2639 0 : ix = incx * (1 - n);
2640 0 : if(incy<0)
2641 0 : iy = incy * (1 - n);
2642 :
2643 0 : for(i=0;i<n;i++,ix+=incx,iy+=incy) {
2644 0 : d1 = dx[ix];
2645 0 : dx[ix] = dy[iy];
2646 0 : dy[iy] = d1;
2647 : }
2648 : }
2649 : return;
2650 : }
2651 :
2652 : }
2653 : }
2654 : #include <cctype>
2655 : #include <cmath>
2656 :
2657 : #include "real.h"
2658 : #include "blas.h"
2659 :
2660 : namespace PLMD{
2661 : namespace blas{
2662 : void
2663 0 : PLUMED_BLAS_F77_FUNC(ssymv,SSYMV)(const char *uplo,
2664 : int *n__,
2665 : float *alpha__,
2666 : float *a,
2667 : int *lda__,
2668 : float *x,
2669 : int *incx__,
2670 : float *beta__,
2671 : float *y,
2672 : int *incy__)
2673 : {
2674 0 : const char ch=std::toupper(*uplo);
2675 : int kx,ky,i,j,ix,iy,jx,jy;
2676 : float temp1,temp2;
2677 :
2678 0 : int n = *n__;
2679 0 : int lda = *lda__;
2680 0 : int incx = *incx__;
2681 0 : int incy = *incy__;
2682 0 : float alpha = *alpha__;
2683 0 : float beta = *beta__;
2684 :
2685 0 : if(n<=0 || incx==0 || incy==0)
2686 : return;
2687 :
2688 0 : if(incx>0)
2689 : kx = 1;
2690 : else
2691 0 : kx = 1 - (n -1)*(incx);
2692 :
2693 0 : if(incy>0)
2694 : ky = 1;
2695 : else
2696 0 : ky = 1 - (n -1)*(incy);
2697 :
2698 0 : if(std::abs(beta-1.0)>PLUMED_GMX_FLOAT_EPS) {
2699 0 : if(incy==1) {
2700 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
2701 0 : for(i=1;i<=n;i++)
2702 0 : y[i-1] = 0.0;
2703 : else
2704 0 : for(i=1;i<=n;i++)
2705 0 : y[i-1] *= beta;
2706 : } else {
2707 : /* non-unit incr. */
2708 : iy = ky;
2709 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
2710 0 : for(i=1;i<=n;i++) {
2711 0 : y[iy-1] = 0.0;
2712 0 : iy += incy;
2713 : }
2714 : else
2715 0 : for(i=1;i<=n;i++) {
2716 0 : y[iy-1] *= beta;
2717 0 : iy += incy;
2718 : }
2719 : }
2720 : }
2721 :
2722 0 : if(std::abs(alpha)<PLUMED_GMX_FLOAT_MIN)
2723 : return;
2724 :
2725 0 : if(ch=='U') {
2726 0 : if(incx==1 && incy==1) {
2727 0 : for(j=1;j<=n;j++) {
2728 0 : temp1 = alpha * x[j-1];
2729 : temp2 = 0.0;
2730 0 : for(i=1;i<j;i++) {
2731 0 : y[i-1] += temp1*a[(j-1)*(lda)+(i-1)];
2732 0 : temp2 += a[(j-1)*(lda)+(i-1)] * x[i-1];
2733 : }
2734 0 : y[j-1] += temp1*a[(j-1)*(lda)+(j-1)] + alpha *temp2;
2735 : }
2736 : } else {
2737 : /* non-unit incr. */
2738 : jx = kx;
2739 : jy = ky;
2740 0 : for(j=1;j<=n;j++) {
2741 0 : temp1 = alpha * x[jx-1];
2742 : temp2 = 0.0;
2743 : ix = kx;
2744 : iy = ky;
2745 0 : for(i=1;i<j;i++) {
2746 0 : y[iy-1] += temp1 * a[(j-1)*(lda)+(i-1)];
2747 0 : temp2 += a[(j-1)*(lda)+(i-1)] * x[ix-1];
2748 0 : ix += incx;
2749 0 : iy += incy;
2750 : }
2751 0 : y[jy-1] += temp1*a[(j-1)*(lda)+(j-1)] + alpha*temp2;
2752 0 : jx += incx;
2753 0 : jy += incy;
2754 : }
2755 : }
2756 : } else {
2757 : /* lower */
2758 0 : if(incx==1 && incy==1) {
2759 0 : for(j=1;j<=n;j++) {
2760 0 : temp1 = alpha * x[j-1];
2761 : temp2 = 0.0;
2762 0 : y[j-1] += temp1 * a[(j-1)*(lda)+(j-1)];
2763 0 : for(i=j+1;i<=n;i++) {
2764 0 : y[i-1] += temp1*a[(j-1)*(lda)+(i-1)];
2765 0 : temp2 += a[(j-1)*(lda)+(i-1)] * x[i-1];
2766 : }
2767 0 : y[j-1] += alpha *temp2;
2768 : }
2769 : } else {
2770 : /* non-unit incr. */
2771 : jx = kx;
2772 : jy = ky;
2773 0 : for(j=1;j<=n;j++) {
2774 0 : temp1 = alpha * x[jx-1];
2775 : temp2 = 0.0;
2776 0 : y[jy-1] += temp1 * a[(j-1)*(lda)+(j-1)];
2777 : ix = jx;
2778 : iy = jy;
2779 0 : for(i=j+1;i<=n;i++) {
2780 0 : ix += incx;
2781 0 : iy += incy;
2782 0 : y[iy-1] += temp1 * a[(j-1)*(lda)+(i-1)];
2783 0 : temp2 += a[(j-1)*(lda)+(i-1)] * x[ix-1];
2784 : }
2785 0 : y[jy-1] += alpha*temp2;
2786 0 : jx += incx;
2787 0 : jy += incy;
2788 : }
2789 : }
2790 : }
2791 : return;
2792 : }
2793 : }
2794 : }
2795 : #include <cctype>
2796 : #include <cmath>
2797 :
2798 : #include "real.h"
2799 : #include "blas.h"
2800 :
2801 : namespace PLMD{
2802 : namespace blas{
2803 : void
2804 0 : PLUMED_BLAS_F77_FUNC(ssyr2,SSYR2)(const char * uplo,
2805 : int * n__,
2806 : float * alpha__,
2807 : float * x,
2808 : int * incx__,
2809 : float * y,
2810 : int * incy__,
2811 : float * a,
2812 : int * lda__)
2813 : {
2814 : int kx,ky,ix,iy,jx,jy,j,i;
2815 : float temp1,temp2;
2816 0 : const char ch=std::toupper(*uplo);
2817 :
2818 0 : int n = *n__;
2819 0 : int lda = *lda__;
2820 0 : int incx = *incx__;
2821 0 : int incy = *incy__;
2822 0 : float alpha = *alpha__;
2823 :
2824 0 : if(n<=0 || std::abs(alpha)<PLUMED_GMX_FLOAT_MIN || incx==0 || incy==0 ||
2825 0 : (ch != 'U' && ch != 'L'))
2826 : return;
2827 :
2828 : jx = jy = kx = ky = 0;
2829 :
2830 : /* init start points for non-unit increments */
2831 0 : if(incx!=1 || incy!=1) {
2832 0 : if(incx>0)
2833 : kx = 1;
2834 : else
2835 0 : kx = 1 - (n - 1)*(incx);
2836 0 : if(incy>0)
2837 : ky = 1;
2838 : else
2839 0 : ky = 1 - (n - 1)*(incy);
2840 :
2841 : jx = kx;
2842 : jy = ky;
2843 : }
2844 :
2845 0 : if(ch == 'U') {
2846 : /* Data in upper part of A */
2847 0 : if(incx==1 && incy==1) {
2848 : /* Unit increments for both x and y */
2849 0 : for(j=1;j<=n;j++) {
2850 0 : if( std::abs(x[j-1])>PLUMED_GMX_FLOAT_MIN || std::abs(y[j-1])>PLUMED_GMX_FLOAT_MIN ) {
2851 0 : temp1 = alpha * y[j-1];
2852 0 : temp2 = alpha * x[j-1];
2853 0 : for(i=1;i<=j;i++)
2854 0 : a[(j-1)*(lda)+(i-1)] += x[i-1]*temp1 + y[i-1]*temp2;
2855 : }
2856 : }
2857 : } else {
2858 :
2859 : /* non-unit increments */
2860 0 : for(j=1;j<=n;j++) {
2861 :
2862 0 : if( std::abs(x[jx-1])>PLUMED_GMX_FLOAT_MIN || std::abs(y[jy-1])>PLUMED_GMX_FLOAT_MIN ) {
2863 0 : temp1 = alpha * y[jy-1];
2864 0 : temp2 = alpha * x[jx-1];
2865 : ix = kx;
2866 : iy = ky;
2867 0 : for(i=1;i<=j;i++) {
2868 0 : a[(j-1)*(lda)+(i-1)] += x[ix-1]*temp1 + y[iy-1]*temp2;
2869 0 : ix += incx;
2870 0 : iy += incy;
2871 : }
2872 : }
2873 0 : jx += incx;
2874 0 : jy += incy;
2875 : }
2876 : }
2877 : } else {
2878 : /* Data in lower part of A */
2879 0 : if(incx==1 && incy==1) {
2880 : /* Unit increments for both x and y */
2881 0 : for(j=1;j<=n;j++) {
2882 0 : if( std::abs(x[j-1])>PLUMED_GMX_FLOAT_MIN || std::abs(y[j-1])>PLUMED_GMX_FLOAT_MIN ) {
2883 0 : temp1 = alpha * y[j-1];
2884 0 : temp2 = alpha * x[j-1];
2885 0 : for(i=j;i<=n;i++)
2886 0 : a[(j-1)*(lda)+(i-1)] += x[i-1]*temp1 + y[i-1]*temp2;
2887 : }
2888 : }
2889 : } else {
2890 :
2891 : /* non-unit increments */
2892 0 : for(j=1;j<=n;j++) {
2893 :
2894 0 : if( std::abs(x[jx-1])>PLUMED_GMX_FLOAT_MIN || std::abs(y[jy-1])>PLUMED_GMX_FLOAT_MIN ) {
2895 0 : temp1 = alpha * y[jy-1];
2896 0 : temp2 = alpha * x[jx-1];
2897 : ix = jx;
2898 : iy = jy;
2899 0 : for(i=j;i<=n;i++) {
2900 0 : a[(j-1)*(lda)+(i-1)] += x[ix-1]*temp1 + y[iy-1]*temp2;
2901 0 : ix += incx;
2902 0 : iy += incy;
2903 : }
2904 : }
2905 0 : jx += incx;
2906 0 : jy += incy;
2907 : }
2908 : }
2909 : }
2910 :
2911 : return;
2912 : }
2913 : }
2914 : }
2915 : #include <cctype>
2916 : #include <cmath>
2917 :
2918 : #include "real.h"
2919 : #include "blas.h"
2920 :
2921 : namespace PLMD{
2922 : namespace blas{
2923 : void
2924 0 : PLUMED_BLAS_F77_FUNC(ssyr2k,SSYR2K)(const char *uplo,
2925 : const char *trans,
2926 : int *n__,
2927 : int *k__,
2928 : float *alpha__,
2929 : float *a,
2930 : int *lda__,
2931 : float *b,
2932 : int *ldb__,
2933 : float *beta__,
2934 : float *c,
2935 : int *ldc__)
2936 : {
2937 : char ch1,ch2;
2938 : int i,j,l;
2939 : float temp1,temp2;
2940 :
2941 0 : int n = *n__;
2942 0 : int k = *k__;
2943 0 : int lda = *lda__;
2944 0 : int ldb = *ldb__;
2945 0 : int ldc = *ldc__;
2946 :
2947 0 : float alpha = *alpha__;
2948 0 : float beta = *beta__;
2949 :
2950 0 : ch1 = std::toupper(*uplo);
2951 0 : ch2 = std::toupper(*trans);
2952 :
2953 0 : if(n==0 || ( ( std::abs(alpha)<PLUMED_GMX_FLOAT_MIN || k==0 ) && std::abs(beta-1.0)<PLUMED_GMX_FLOAT_EPS))
2954 : return;
2955 :
2956 0 : if(std::abs(alpha)<PLUMED_GMX_FLOAT_MIN) {
2957 0 : if(ch1=='U') {
2958 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
2959 0 : for(j=1;j<=n;j++)
2960 0 : for(i=1;i<=j;i++)
2961 0 : c[(j-1)*(ldc)+(i-1)] = 0.0;
2962 : else
2963 0 : for(j=1;j<=n;j++)
2964 0 : for(i=1;i<=j;i++)
2965 0 : c[(j-1)*(ldc)+(i-1)] *= beta;
2966 : } else {
2967 : /* lower */
2968 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
2969 0 : for(j=1;j<=n;j++)
2970 0 : for(i=j;i<=n;i++)
2971 0 : c[(j-1)*(ldc)+(i-1)] = 0.0;
2972 : else
2973 0 : for(j=1;j<=n;j++)
2974 0 : for(i=j;i<=n;i++)
2975 0 : c[(j-1)*(ldc)+(i-1)] *= beta;
2976 : }
2977 0 : return;
2978 : }
2979 :
2980 0 : if(ch2=='N') {
2981 0 : if(ch1=='U') {
2982 0 : for(j=1;j<=n;j++) {
2983 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
2984 0 : for(i=1;i<=j;i++)
2985 0 : c[(j-1)*(ldc)+(i-1)] = 0.0;
2986 0 : else if(std::abs(beta-1.0)>PLUMED_GMX_FLOAT_EPS)
2987 0 : for(i=1;i<=j;i++)
2988 0 : c[(j-1)*(ldc)+(i-1)] *= beta;
2989 0 : for(l=1;l<=k;l++) {
2990 0 : if( std::abs(a[(l-1)*(lda)+(j-1)])>PLUMED_GMX_FLOAT_MIN ||
2991 0 : std::abs(b[(l-1)*(ldb)+(j-1)])>PLUMED_GMX_FLOAT_MIN) {
2992 0 : temp1 = alpha * b[(l-1)*(ldb)+(j-1)];
2993 0 : temp2 = alpha * a[(l-1)*(lda)+(j-1)];
2994 0 : for(i=1;i<=j;i++)
2995 0 : c[(j-1)*(ldc)+(i-1)] +=
2996 0 : a[(l-1)*(lda)+(i-1)] * temp1 +
2997 0 : b[(l-1)*(ldb)+(i-1)] * temp2;
2998 : }
2999 : }
3000 : }
3001 : } else {
3002 : /* lower */
3003 0 : for(j=1;j<=n;j++) {
3004 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
3005 0 : for(i=j;i<=n;i++)
3006 0 : c[(j-1)*(ldc)+(i-1)] = 0.0;
3007 0 : else if(std::abs(beta-1.0)>PLUMED_GMX_FLOAT_EPS)
3008 0 : for(i=j;i<=n;i++)
3009 0 : c[(j-1)*(ldc)+(i-1)] *= beta;
3010 0 : for(l=1;l<=k;l++) {
3011 0 : if( std::abs(a[(l-1)*(lda)+(j-1)])>PLUMED_GMX_FLOAT_MIN ||
3012 0 : std::abs(b[(l-1)*(ldb)+(j-1)])>PLUMED_GMX_FLOAT_MIN) {
3013 0 : temp1 = alpha * b[(l-1)*(ldb)+(j-1)];
3014 0 : temp2 = alpha * a[(l-1)*(lda)+(j-1)];
3015 0 : for(i=j;i<=n;i++)
3016 0 : c[(j-1)*(ldc)+(i-1)] +=
3017 0 : a[(l-1)*(lda)+(i-1)] * temp1 +
3018 0 : b[(l-1)*(ldb)+(i-1)] * temp2;
3019 : }
3020 : }
3021 : }
3022 : }
3023 : } else {
3024 : /* transpose */
3025 0 : if(ch1=='U') {
3026 0 : for(j=1;j<=n;j++)
3027 0 : for(i=1;i<=j;i++) {
3028 : temp1 = 0.0;
3029 : temp2 = 0.0;
3030 0 : for (l=1;l<=k;l++) {
3031 0 : temp1 += a[(i-1)*(lda)+(l-1)] * b[(j-1)*(ldb)+(l-1)];
3032 0 : temp2 += b[(i-1)*(ldb)+(l-1)] * a[(j-1)*(lda)+(l-1)];
3033 : }
3034 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
3035 0 : c[(j-1)*(ldc)+(i-1)] = alpha * (temp1 + temp2);
3036 : else
3037 0 : c[(j-1)*(ldc)+(i-1)] = beta * c[(j-1)*(ldc)+(i-1)] +
3038 0 : alpha * (temp1 + temp2);
3039 : }
3040 : } else {
3041 : /* lower */
3042 0 : for(j=1;j<=n;j++)
3043 0 : for(i=j;i<=n;i++) {
3044 : temp1 = 0.0;
3045 : temp2 = 0.0;
3046 0 : for (l=1;l<=k;l++) {
3047 0 : temp1 += a[(i-1)*(lda)+(l-1)] * b[(j-1)*(ldb)+(l-1)];
3048 0 : temp2 += b[(i-1)*(ldb)+(l-1)] * a[(j-1)*(lda)+(l-1)];
3049 : }
3050 0 : if(std::abs(beta)<PLUMED_GMX_FLOAT_MIN)
3051 0 : c[(j-1)*(ldc)+(i-1)] = alpha * (temp1 + temp2);
3052 : else
3053 0 : c[(j-1)*(ldc)+(i-1)] = beta * c[(j-1)*(ldc)+(i-1)] +
3054 0 : alpha * (temp1 + temp2);
3055 : }
3056 : }
3057 : }
3058 : return;
3059 : }
3060 : }
3061 : }
3062 : #include <cmath>
3063 :
3064 : #include "real.h"
3065 : #include "blas.h"
3066 :
3067 : namespace PLMD{
3068 : namespace blas{
3069 : void
3070 0 : PLUMED_BLAS_F77_FUNC(strmm,STRMM)(const char *side,
3071 : const char *uplo,
3072 : const char *transa,
3073 : const char *diag,
3074 : int *m__,
3075 : int *n__,
3076 : float *alpha__,
3077 : float *a,
3078 : int *lda__,
3079 : float *b,
3080 : int *ldb__)
3081 : {
3082 : int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
3083 :
3084 0 : int m = *m__;
3085 0 : int n = *n__;
3086 0 : int lda = *lda__;
3087 0 : int ldb = *ldb__;
3088 0 : float alpha = *alpha__;
3089 :
3090 : /* Local variables */
3091 : int i__, j, k;
3092 : float temp;
3093 : int lside;
3094 : int upper;
3095 : int nounit;
3096 : a_dim1 = lda;
3097 0 : a_offset = 1 + a_dim1;
3098 0 : a -= a_offset;
3099 : b_dim1 = ldb;
3100 0 : b_offset = 1 + b_dim1;
3101 0 : b -= b_offset;
3102 :
3103 : /* Function Body */
3104 0 : lside = (*side=='L' || *side=='l');
3105 :
3106 0 : nounit = (*diag=='N' || *diag=='n');
3107 0 : upper = (*uplo=='U' || *uplo=='u');
3108 :
3109 0 : if (n == 0) {
3110 : return;
3111 : }
3112 0 : if (std::abs(alpha)<PLUMED_GMX_FLOAT_MIN) {
3113 : i__1 = n;
3114 0 : for (j = 1; j <= i__1; ++j) {
3115 : i__2 = m;
3116 0 : for (i__ = 1; i__ <= i__2; ++i__) {
3117 0 : b[i__ + j * b_dim1] = 0.;
3118 : }
3119 : }
3120 : return;
3121 : }
3122 0 : if (lside) {
3123 0 : if (*transa=='N' || *transa=='n') {
3124 0 : if (upper) {
3125 : i__1 = n;
3126 0 : for (j = 1; j <= i__1; ++j) {
3127 : i__2 = m;
3128 0 : for (k = 1; k <= i__2; ++k) {
3129 0 : if ( std::abs(b[k + j * b_dim1])>PLUMED_GMX_FLOAT_MIN) {
3130 0 : temp = alpha * b[k + j * b_dim1];
3131 : i__3 = k - 1;
3132 0 : for (i__ = 1; i__ <= i__3; ++i__) {
3133 0 : b[i__ + j * b_dim1] += temp * a[i__ + k *
3134 0 : a_dim1];
3135 : }
3136 0 : if (nounit) {
3137 0 : temp *= a[k + k * a_dim1];
3138 : }
3139 0 : b[k + j * b_dim1] = temp;
3140 : }
3141 : }
3142 : }
3143 : } else {
3144 : i__1 = n;
3145 0 : for (j = 1; j <= i__1; ++j) {
3146 0 : for (k = m; k >= 1; --k) {
3147 0 : if (std::abs(b[k + j * b_dim1])>PLUMED_GMX_FLOAT_MIN) {
3148 0 : temp = alpha * b[k + j * b_dim1];
3149 0 : b[k + j * b_dim1] = temp;
3150 0 : if (nounit) {
3151 0 : b[k + j * b_dim1] *= a[k + k * a_dim1];
3152 : }
3153 : i__2 = m;
3154 0 : for (i__ = k + 1; i__ <= i__2; ++i__) {
3155 0 : b[i__ + j * b_dim1] += temp * a[i__ + k *
3156 0 : a_dim1];
3157 : }
3158 : }
3159 : }
3160 : }
3161 : }
3162 : } else {
3163 :
3164 0 : if (upper) {
3165 : i__1 = n;
3166 0 : for (j = 1; j <= i__1; ++j) {
3167 0 : for (i__ = m; i__ >= 1; --i__) {
3168 0 : temp = b[i__ + j * b_dim1];
3169 0 : if (nounit) {
3170 0 : temp *= a[i__ + i__ * a_dim1];
3171 : }
3172 : i__2 = i__ - 1;
3173 0 : for (k = 1; k <= i__2; ++k) {
3174 0 : temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
3175 : }
3176 0 : b[i__ + j * b_dim1] = alpha * temp;
3177 : }
3178 : }
3179 : } else {
3180 : i__1 = n;
3181 0 : for (j = 1; j <= i__1; ++j) {
3182 : i__2 = m;
3183 0 : for (i__ = 1; i__ <= i__2; ++i__) {
3184 0 : temp = b[i__ + j * b_dim1];
3185 0 : if (nounit) {
3186 0 : temp *= a[i__ + i__ * a_dim1];
3187 : }
3188 : i__3 = m;
3189 0 : for (k = i__ + 1; k <= i__3; ++k) {
3190 0 : temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
3191 : }
3192 0 : b[i__ + j * b_dim1] = alpha * temp;
3193 : }
3194 : }
3195 : }
3196 : }
3197 : } else {
3198 0 : if (*transa=='N' || *transa=='n') {
3199 :
3200 0 : if (upper) {
3201 0 : for (j = n; j >= 1; --j) {
3202 : temp = alpha;
3203 0 : if (nounit) {
3204 0 : temp *= a[j + j * a_dim1];
3205 : }
3206 : i__1 = m;
3207 0 : for (i__ = 1; i__ <= i__1; ++i__) {
3208 0 : b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
3209 : }
3210 : i__1 = j - 1;
3211 0 : for (k = 1; k <= i__1; ++k) {
3212 0 : if ( std::abs(a[k + j * a_dim1])>PLUMED_GMX_FLOAT_MIN) {
3213 0 : temp = alpha * a[k + j * a_dim1];
3214 : i__2 = m;
3215 0 : for (i__ = 1; i__ <= i__2; ++i__) {
3216 0 : b[i__ + j * b_dim1] += temp * b[i__ + k *
3217 0 : b_dim1];
3218 : }
3219 : }
3220 : }
3221 : }
3222 : } else {
3223 : i__1 = n;
3224 0 : for (j = 1; j <= i__1; ++j) {
3225 : temp = alpha;
3226 0 : if (nounit) {
3227 0 : temp *= a[j + j * a_dim1];
3228 : }
3229 : i__2 = m;
3230 0 : for (i__ = 1; i__ <= i__2; ++i__) {
3231 0 : b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
3232 : }
3233 : i__2 = n;
3234 0 : for (k = j + 1; k <= i__2; ++k) {
3235 0 : if ( std::abs(a[k + j * a_dim1])>PLUMED_GMX_FLOAT_MIN) {
3236 0 : temp = alpha * a[k + j * a_dim1];
3237 : i__3 = m;
3238 0 : for (i__ = 1; i__ <= i__3; ++i__) {
3239 0 : b[i__ + j * b_dim1] += temp * b[i__ + k *
3240 0 : b_dim1];
3241 : }
3242 : }
3243 : }
3244 : }
3245 : }
3246 : } else {
3247 :
3248 0 : if (upper) {
3249 : i__1 = n;
3250 0 : for (k = 1; k <= i__1; ++k) {
3251 : i__2 = k - 1;
3252 0 : for (j = 1; j <= i__2; ++j) {
3253 0 : if ( std::abs(a[j + k * a_dim1])>PLUMED_GMX_FLOAT_MIN ) {
3254 0 : temp = alpha * a[j + k * a_dim1];
3255 : i__3 = m;
3256 0 : for (i__ = 1; i__ <= i__3; ++i__) {
3257 0 : b[i__ + j * b_dim1] += temp * b[i__ + k *
3258 0 : b_dim1];
3259 : }
3260 : }
3261 : }
3262 : temp = alpha;
3263 0 : if (nounit) {
3264 0 : temp *= a[k + k * a_dim1];
3265 : }
3266 0 : if ( std::abs(temp-1.0)>PLUMED_GMX_FLOAT_EPS) {
3267 : i__2 = m;
3268 0 : for (i__ = 1; i__ <= i__2; ++i__) {
3269 0 : b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
3270 : }
3271 : }
3272 : }
3273 : } else {
3274 0 : for (k = n; k >= 1; --k) {
3275 : i__1 = n;
3276 0 : for (j = k + 1; j <= i__1; ++j) {
3277 0 : if ( std::abs(a[j + k * a_dim1])>PLUMED_GMX_FLOAT_MIN) {
3278 0 : temp = alpha * a[j + k * a_dim1];
3279 : i__2 = m;
3280 0 : for (i__ = 1; i__ <= i__2; ++i__) {
3281 0 : b[i__ + j * b_dim1] += temp * b[i__ + k *
3282 0 : b_dim1];
3283 : }
3284 : }
3285 : }
3286 : temp = alpha;
3287 0 : if (nounit) {
3288 0 : temp *= a[k + k * a_dim1];
3289 : }
3290 0 : if ( std::abs(temp-1.0)>PLUMED_GMX_FLOAT_EPS) {
3291 : i__1 = m;
3292 0 : for (i__ = 1; i__ <= i__1; ++i__) {
3293 0 : b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
3294 : }
3295 : }
3296 : }
3297 : }
3298 : }
3299 : }
3300 :
3301 : return;
3302 :
3303 : }
3304 :
3305 :
3306 : }
3307 : }
3308 : #include <cmath>
3309 :
3310 : #include "real.h"
3311 : #include "blas.h"
3312 :
3313 : namespace PLMD{
3314 : namespace blas{
3315 : void
3316 0 : PLUMED_BLAS_F77_FUNC(strmv,STRMV)(const char *uplo,
3317 : const char *trans,
3318 : const char *diag,
3319 : int *n__,
3320 : float *a,
3321 : int *lda__,
3322 : float *x,
3323 : int *incx__)
3324 : {
3325 : int a_dim1, a_offset, i__1, i__2;
3326 :
3327 : int i__, j, ix, jx, kx;
3328 : float temp;
3329 : int nounit;
3330 :
3331 0 : int n = *n__;
3332 0 : int lda = *lda__;
3333 0 : int incx = *incx__;
3334 :
3335 : a_dim1 = lda;
3336 0 : a_offset = 1 + a_dim1;
3337 0 : a -= a_offset;
3338 0 : --x;
3339 :
3340 0 : if (n == 0) {
3341 : return;
3342 : }
3343 :
3344 0 : nounit = (*diag=='n' || *diag=='N');
3345 :
3346 0 : if (incx <= 0) {
3347 0 : kx = 1 - (n - 1) * incx;
3348 : } else {
3349 : kx = 1;
3350 : }
3351 :
3352 0 : if (*trans=='N' || *trans=='n') {
3353 :
3354 0 : if (*uplo=='U' || *uplo=='u') {
3355 0 : if (incx == 1) {
3356 : i__1 = n;
3357 0 : for (j = 1; j <= i__1; ++j) {
3358 0 : if (std::abs(x[j])>PLUMED_GMX_FLOAT_MIN) {
3359 : temp = x[j];
3360 : i__2 = j - 1;
3361 0 : for (i__ = 1; i__ <= i__2; ++i__) {
3362 0 : x[i__] += temp * a[i__ + j * a_dim1];
3363 : }
3364 0 : if (nounit) {
3365 0 : x[j] *= a[j + j * a_dim1];
3366 : }
3367 : }
3368 : }
3369 : } else {
3370 : jx = kx;
3371 : i__1 = n;
3372 0 : for (j = 1; j <= i__1; ++j) {
3373 0 : if (std::abs(x[jx])>PLUMED_GMX_FLOAT_MIN) {
3374 : temp = x[jx];
3375 : ix = kx;
3376 : i__2 = j - 1;
3377 0 : for (i__ = 1; i__ <= i__2; ++i__) {
3378 0 : x[ix] += temp * a[i__ + j * a_dim1];
3379 0 : ix += incx;
3380 : }
3381 0 : if (nounit) {
3382 0 : x[jx] *= a[j + j * a_dim1];
3383 : }
3384 : }
3385 0 : jx += incx;
3386 : }
3387 : }
3388 : } else {
3389 0 : if (incx == 1) {
3390 0 : for (j = n; j >= 1; --j) {
3391 0 : if (std::abs(x[j])>PLUMED_GMX_FLOAT_MIN) {
3392 : temp = x[j];
3393 : i__1 = j + 1;
3394 0 : for (i__ = n; i__ >= i__1; --i__) {
3395 0 : x[i__] += temp * a[i__ + j * a_dim1];
3396 : }
3397 0 : if (nounit) {
3398 0 : x[j] *= a[j + j * a_dim1];
3399 : }
3400 : }
3401 : }
3402 : } else {
3403 0 : kx += (n - 1) * incx;
3404 : jx = kx;
3405 0 : for (j = n; j >= 1; --j) {
3406 0 : if (std::abs(x[jx])>PLUMED_GMX_FLOAT_MIN) {
3407 : temp = x[jx];
3408 : ix = kx;
3409 : i__1 = j + 1;
3410 0 : for (i__ = n; i__ >= i__1; --i__) {
3411 0 : x[ix] += temp * a[i__ + j * a_dim1];
3412 0 : ix -= incx;
3413 : }
3414 0 : if (nounit) {
3415 0 : x[jx] *= a[j + j * a_dim1];
3416 : }
3417 : }
3418 0 : jx -= incx;
3419 : }
3420 : }
3421 : }
3422 : } else {
3423 :
3424 0 : if (*uplo=='U' || *uplo=='u') {
3425 0 : if (incx == 1) {
3426 0 : for (j = n; j >= 1; --j) {
3427 0 : temp = x[j];
3428 0 : if (nounit) {
3429 0 : temp *= a[j + j * a_dim1];
3430 : }
3431 0 : for (i__ = j - 1; i__ >= 1; --i__) {
3432 0 : temp += a[i__ + j * a_dim1] * x[i__];
3433 : }
3434 0 : x[j] = temp;
3435 : }
3436 : } else {
3437 0 : jx = kx + (n - 1) * incx;
3438 0 : for (j = n; j >= 1; --j) {
3439 0 : temp = x[jx];
3440 : ix = jx;
3441 0 : if (nounit) {
3442 0 : temp *= a[j + j * a_dim1];
3443 : }
3444 0 : for (i__ = j - 1; i__ >= 1; --i__) {
3445 0 : ix -= incx;
3446 0 : temp += a[i__ + j * a_dim1] * x[ix];
3447 : }
3448 0 : x[jx] = temp;
3449 0 : jx -= incx;
3450 : }
3451 : }
3452 : } else {
3453 0 : if (incx == 1) {
3454 : i__1 = n;
3455 0 : for (j = 1; j <= i__1; ++j) {
3456 0 : temp = x[j];
3457 0 : if (nounit) {
3458 0 : temp *= a[j + j * a_dim1];
3459 : }
3460 : i__2 = n;
3461 0 : for (i__ = j + 1; i__ <= i__2; ++i__) {
3462 0 : temp += a[i__ + j * a_dim1] * x[i__];
3463 : }
3464 0 : x[j] = temp;
3465 : }
3466 : } else {
3467 : jx = kx;
3468 : i__1 = n;
3469 0 : for (j = 1; j <= i__1; ++j) {
3470 0 : temp = x[jx];
3471 : ix = jx;
3472 0 : if (nounit) {
3473 0 : temp *= a[j + j * a_dim1];
3474 : }
3475 : i__2 = n;
3476 0 : for (i__ = j + 1; i__ <= i__2; ++i__) {
3477 0 : ix += incx;
3478 0 : temp += a[i__ + j * a_dim1] * x[ix];
3479 : }
3480 0 : x[jx] = temp;
3481 0 : jx += incx;
3482 : }
3483 : }
3484 : }
3485 : }
3486 :
3487 : return;
3488 :
3489 : }
3490 :
3491 :
3492 : }
3493 : }
3494 : #include <cctype>
3495 : #include <cmath>
3496 :
3497 : #include "real.h"
3498 : #include "blas.h"
3499 :
3500 : namespace PLMD{
3501 : namespace blas{
3502 : void
3503 0 : PLUMED_BLAS_F77_FUNC(strsm,STRSM)(const char * side,
3504 : const char * uplo,
3505 : const char * transa,
3506 : const char * diag,
3507 : int * m__,
3508 : int * n__,
3509 : float *alpha__,
3510 : float *a,
3511 : int * lda__,
3512 : float *b,
3513 : int * ldb__)
3514 : {
3515 0 : const char xside = std::toupper(*side);
3516 0 : const char xuplo = std::toupper(*uplo);
3517 0 : const char xtrans = std::toupper(*transa);
3518 0 : const char xdiag = std::toupper(*diag);
3519 : int i,j,k;
3520 : float temp;
3521 :
3522 0 : int m = *m__;
3523 0 : int n = *n__;
3524 0 : int lda = *lda__;
3525 0 : int ldb = *ldb__;
3526 0 : float alpha = *alpha__;
3527 :
3528 0 : if(n<=0)
3529 : return;
3530 :
3531 :
3532 0 : if(std::abs(alpha)<PLUMED_GMX_FLOAT_MIN) {
3533 0 : for(j=0;j<n;j++)
3534 0 : for(i=0;i<m;i++)
3535 0 : b[j*(ldb)+i] = 0.0;
3536 : return;
3537 : }
3538 :
3539 0 : if(xside=='L') {
3540 : /* left side */
3541 0 : if(xtrans=='N') {
3542 : /* No transpose */
3543 0 : if(xuplo=='U') {
3544 : /* upper */
3545 0 : for(j=0;j<n;j++) {
3546 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_FLOAT_EPS) {
3547 0 : for(i=0;i<m;i++)
3548 0 : b[j*(ldb)+i] *= alpha;
3549 : }
3550 0 : for(k=m-1;k>=0;k--) {
3551 0 : if( std::abs(b[j*(ldb)+k])>PLUMED_GMX_FLOAT_MIN) {
3552 0 : if(xdiag=='N')
3553 0 : b[j*(ldb)+k] /= a[k*(lda)+k];
3554 0 : for(i=0;i<k;i++)
3555 0 : b[j*(ldb)+i] -= b[j*(ldb)+k]*a[k*(lda)+i];
3556 : }
3557 : }
3558 : }
3559 : } else {
3560 : /* lower */
3561 0 : for(j=0;j<n;j++) {
3562 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_FLOAT_EPS)
3563 0 : for(i=0;i<m;i++)
3564 0 : b[j*(ldb)+i] *= alpha;
3565 0 : for(k=0;k<m;k++) {
3566 0 : if( std::abs(b[j*(ldb)+k])>PLUMED_GMX_FLOAT_MIN) {
3567 0 : if(xdiag=='N')
3568 0 : b[j*(ldb)+k] /= a[k*(lda)+k];
3569 0 : for(i=k+1;i<m;i++)
3570 0 : b[j*(ldb)+i] -= b[j*(ldb)+k]*a[k*(lda)+i];
3571 : }
3572 : }
3573 : }
3574 : }
3575 : } else {
3576 : /* Transpose */
3577 0 : if(xuplo=='U') {
3578 : /* upper */
3579 0 : for(j=0;j<n;j++) {
3580 0 : for(i=0;i<m;i++) {
3581 0 : temp = alpha * b[j*(ldb)+i];
3582 0 : for(k=0;k<i;k++)
3583 0 : temp -= a[i*(lda)+k] * b[j*(ldb)+k];
3584 0 : if(xdiag=='N')
3585 0 : temp /= a[i*(lda)+i];
3586 0 : b[j*(ldb)+i] = temp;
3587 : }
3588 : }
3589 : } else {
3590 : /* lower */
3591 0 : for(j=0;j<n;j++) {
3592 0 : for(i=m-1;i>=0;i--) {
3593 0 : temp = alpha * b[j*(ldb)+i];
3594 0 : for(k=i+1;k<m;k++)
3595 0 : temp -= a[i*(lda)+k] * b[j*(ldb)+k];
3596 0 : if(xdiag=='N')
3597 0 : temp /= a[i*(lda)+i];
3598 0 : b[j*(ldb)+i] = temp;
3599 : }
3600 : }
3601 : }
3602 : }
3603 : } else {
3604 : /* right side */
3605 0 : if(xtrans=='N') {
3606 : /* No transpose */
3607 0 : if(xuplo=='U') {
3608 : /* upper */
3609 0 : for(j=0;j<n;j++) {
3610 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_FLOAT_EPS)
3611 0 : for(i=0;i<m;i++)
3612 0 : b[j*(ldb)+i] *= alpha;
3613 0 : for(k=0;k<j;k++) {
3614 0 : if( std::abs(a[j*(lda)+k])>PLUMED_GMX_FLOAT_MIN) {
3615 0 : for(i=0;i<m;i++)
3616 0 : b[j*(ldb)+i] -= a[j*(lda)+k]*b[k*(ldb)+i];
3617 : }
3618 : }
3619 0 : if(xdiag=='N') {
3620 0 : temp = 1.0/a[j*(lda)+j];
3621 0 : for(i=0;i<m;i++)
3622 0 : b[j*(ldb)+i] *= temp;
3623 : }
3624 : }
3625 : } else {
3626 : /* lower */
3627 0 : for(j=n-1;j>=0;j--) {
3628 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_FLOAT_EPS)
3629 0 : for(i=0;i<m;i++)
3630 0 : b[j*(ldb)+i] *= alpha;
3631 0 : for(k=j+1;k<n;k++) {
3632 0 : if( std::abs(a[j*(lda)+k])>PLUMED_GMX_FLOAT_MIN ) {
3633 0 : for(i=0;i<m;i++)
3634 0 : b[j*(ldb)+i] -= a[j*(lda)+k]*b[k*(ldb)+i];
3635 : }
3636 : }
3637 0 : if(xdiag=='N') {
3638 0 : temp = 1.0/a[j*(lda)+j];
3639 0 : for(i=0;i<m;i++)
3640 0 : b[j*(ldb)+i] *= temp;
3641 : }
3642 : }
3643 : }
3644 : } else {
3645 : /* Transpose */
3646 0 : if(xuplo=='U') {
3647 : /* upper */
3648 0 : for(k=n-1;k>=0;k--) {
3649 0 : if(xdiag=='N') {
3650 0 : temp = 1.0/a[k*(lda)+k];
3651 0 : for(i=0;i<m;i++)
3652 0 : b[k*(ldb)+i] *= temp;
3653 : }
3654 0 : for(j=0;j<k;j++) {
3655 0 : if( std::abs(a[k*(lda)+j])>PLUMED_GMX_FLOAT_MIN) {
3656 : temp = a[k*(lda)+j];
3657 0 : for(i=0;i<m;i++)
3658 0 : b[j*(ldb)+i] -= temp * b[k*(ldb)+i];
3659 : }
3660 : }
3661 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_FLOAT_EPS)
3662 0 : for(i=0;i<m;i++)
3663 0 : b[k*(ldb)+i] *= alpha;
3664 : }
3665 : } else {
3666 : /* lower */
3667 0 : for(k=0;k<n;k++) {
3668 0 : if(xdiag=='N') {
3669 0 : temp = 1.0/a[k*(lda)+k];
3670 0 : for(i=0;i<m;i++)
3671 0 : b[k*(ldb)+i] *= temp;
3672 : }
3673 0 : for(j=k+1;j<n;j++) {
3674 0 : if( std::abs(a[k*(lda)+j])>PLUMED_GMX_FLOAT_MIN) {
3675 : temp = a[k*(lda)+j];
3676 0 : for(i=0;i<m;i++)
3677 0 : b[j*(ldb)+i] -= temp * b[k*(ldb)+i];
3678 : }
3679 : }
3680 0 : if(std::abs(alpha-1.0)>PLUMED_GMX_FLOAT_EPS)
3681 0 : for(i=0;i<m;i++)
3682 0 : b[k*(ldb)+i] *= alpha;
3683 : }
3684 : }
3685 : }
3686 : }
3687 : }
3688 : }
3689 : }
3690 : #endif
|
