ifa: Fixup
[openmx:openmx.git] / src / omxBLAS.f
1
2
3       SUBROUTINE OMXUNSAFEDGEMM (TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,
4      $                   BETA, C, LDC )
5 *     .. Scalar Arguments ..
6       CHARACTER*1        TRANSA, TRANSB
7       INTEGER            M, N, K, LDA, LDB, LDC
8       DOUBLE PRECISION   ALPHA, BETA
9 *     .. Array Arguments ..
10       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * )
11 *     ..
12 *
13 *  Purpose
14 *  =======
15 *
16 *  DGEMM  performs one of the matrix-matrix operations
17 *
18 *     C := alpha*op( A )*op( B ) + beta*C,
19 *
20 *  where  op( X ) is one of
21 *
22 *     op( X ) = X   or   op( X ) = X',
23 *
24 *  alpha and beta are scalars, and A, B and C are matrices, with op( A )
25 *  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
26 *
27 *  Parameters
28 *  ==========
29 *
30 *  TRANSA - CHARACTER*1.
31 *           On entry, TRANSA specifies the form of op( A ) to be used in
32 *           the matrix multiplication as follows:
33 *
34 *              TRANSA = 'N' or 'n',  op( A ) = A.
35 *
36 *              TRANSA = 'T' or 't',  op( A ) = A'.
37 *
38 *              TRANSA = 'C' or 'c',  op( A ) = A'.
39 *
40 *           Unchanged on exit.
41 *
42 *  TRANSB - CHARACTER*1.
43 *           On entry, TRANSB specifies the form of op( B ) to be used in
44 *           the matrix multiplication as follows:
45 *
46 *              TRANSB = 'N' or 'n',  op( B ) = B.
47 *
48 *              TRANSB = 'T' or 't',  op( B ) = B'.
49 *
50 *              TRANSB = 'C' or 'c',  op( B ) = B'.
51 *
52 *           Unchanged on exit.
53 *
54 *  M      - INTEGER.
55 *           On entry,  M  specifies  the number  of rows  of the  matrix
56 *           op( A )  and of the  matrix  C.  M  must  be at least  zero.
57 *           Unchanged on exit.
58 *
59 *  N      - INTEGER.
60 *           On entry,  N  specifies the number  of columns of the matrix
61 *           op( B ) and the number of columns of the matrix C. N must be
62 *           at least zero.
63 *           Unchanged on exit.
64 *
65 *  K      - INTEGER.
66 *           On entry,  K  specifies  the number of columns of the matrix
67 *           op( A ) and the number of rows of the matrix op( B ). K must
68 *           be at least  zero.
69 *           Unchanged on exit.
70 *
71 *  ALPHA  - DOUBLE PRECISION.
72 *           On entry, ALPHA specifies the scalar alpha.
73 *           Unchanged on exit.
74 *
75 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
76 *           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
77 *           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
78 *           part of the array  A  must contain the matrix  A,  otherwise
79 *           the leading  k by m  part of the array  A  must contain  the
80 *           matrix A.
81 *           Unchanged on exit.
82 *
83 *  LDA    - INTEGER.
84 *           On entry, LDA specifies the first dimension of A as declared
85 *           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
86 *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
87 *           least  max( 1, k ).
88 *           Unchanged on exit.
89 *
90 *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
91 *           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
92 *           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
93 *           part of the array  B  must contain the matrix  B,  otherwise
94 *           the leading  n by k  part of the array  B  must contain  the
95 *           matrix B.
96 *           Unchanged on exit.
97 *
98 *  LDB    - INTEGER.
99 *           On entry, LDB specifies the first dimension of B as declared
100 *           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
101 *           LDB must be at least  max( 1, k ), otherwise  LDB must be at
102 *           least  max( 1, n ).
103 *           Unchanged on exit.
104 *
105 *  BETA   - DOUBLE PRECISION.
106 *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
107 *           supplied as zero then C need not be set on input.
108 *           Unchanged on exit.
109 *
110 *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
111 *           Before entry, the leading  m by n  part of the array  C must
112 *           contain the matrix  C,  except when  beta  is zero, in which
113 *           case C need not be set on entry.
114 *           On exit, the array  C  is overwritten by the  m by n  matrix
115 *           ( alpha*op( A )*op( B ) + beta*C ).
116 *
117 *  LDC    - INTEGER.
118 *           On entry, LDC specifies the first dimension of C as declared
119 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
120 *           max( 1, m ).
121 *           Unchanged on exit.
122 *
123 *
124 *  Level 3 Blas routine.
125 *
126 *  -- Written on 8-February-1989.
127 *     Jack Dongarra, Argonne National Laboratory.
128 *     Iain Duff, AERE Harwell.
129 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
130 *     Sven Hammarling, Numerical Algorithms Group Ltd.
131 *
132 *
133 *     .. External Functions ..
134       LOGICAL            LSAME
135       EXTERNAL           LSAME
136 *     .. External Subroutines ..
137       EXTERNAL           XERBLA
138 *     .. Intrinsic Functions ..
139       INTRINSIC          MAX
140 *     .. Local Scalars ..
141       LOGICAL            NOTA, NOTB
142       INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
143       DOUBLE PRECISION   TEMP
144 *     .. Parameters ..
145       DOUBLE PRECISION   ONE         , ZERO
146       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
147 *     ..
148 *     .. Executable Statements ..
149 *
150 *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
151 *     transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
152 *     and  columns of  A  and the  number of  rows  of  B  respectively.
153 *
154       NOTA  = LSAME( TRANSA, 'N' )
155       NOTB  = LSAME( TRANSB, 'N' )
156       IF( NOTA )THEN
157          NROWA = M
158          NCOLA = K
159       ELSE
160          NROWA = K
161          NCOLA = M
162       END IF
163       IF( NOTB )THEN
164          NROWB = K
165       ELSE
166          NROWB = N
167       END IF
168 *
169 *     Test the input parameters.
170 *
171       INFO = 0
172       IF(      ( .NOT.NOTA                 ).AND.
173      $         ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
174      $         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
175          INFO = 1
176       ELSE IF( ( .NOT.NOTB                 ).AND.
177      $         ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
178      $         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
179          INFO = 2
180       ELSE IF( M  .LT.0               )THEN
181          INFO = 3
182       ELSE IF( N  .LT.0               )THEN
183          INFO = 4
184       ELSE IF( K  .LT.0               )THEN
185          INFO = 5
186       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
187          INFO = 8
188       ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
189          INFO = 10
190       ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
191          INFO = 13
192       END IF
193       IF( INFO.NE.0 )THEN
194          CALL XERBLA( 'OMXUNSAFEDGEMM ', INFO )
195          RETURN
196       END IF
197 *
198 *     Quick return if possible.
199 *
200       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
201      $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
202      $   RETURN
203 *
204 *     And if  alpha.eq.zero.
205 *
206       IF( ALPHA.EQ.ZERO )THEN
207          IF( BETA.EQ.ZERO )THEN
208             DO 20, J = 1, N
209                DO 10, I = 1, M
210                   C( I, J ) = ZERO
211    10          CONTINUE
212    20       CONTINUE
213          ELSE
214             DO 40, J = 1, N
215                DO 30, I = 1, M
216                   C( I, J ) = BETA*C( I, J )
217    30          CONTINUE
218    40       CONTINUE
219          END IF
220          RETURN
221       END IF
222 *
223 *     Start the operations.
224 *
225       IF( NOTB )THEN
226          IF( NOTA )THEN
227 *
228 *           Form  C := alpha*A*B + beta*C.
229 *
230             DO 90, J = 1, N
231                IF( BETA.EQ.ZERO )THEN
232                   DO 50, I = 1, M
233                      C( I, J ) = ZERO
234    50             CONTINUE
235                ELSE IF( BETA.NE.ONE )THEN
236                   DO 60, I = 1, M
237                      C( I, J ) = BETA*C( I, J )
238    60             CONTINUE
239                END IF
240                DO 80, L = 1, K
241                   IF( B( L, J ).NE.ZERO )THEN
242                      TEMP = ALPHA*B( L, J )
243                      DO 70, I = 1, M
244                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
245    70                CONTINUE
246                   END IF
247    80          CONTINUE
248    90       CONTINUE
249          ELSE
250 *
251 *           Form  C := alpha*A'*B + beta*C
252 *
253             DO 120, J = 1, N
254                DO 110, I = 1, M
255                   TEMP = ZERO
256                  DO 100, L = 1, K
257                      TEMP = TEMP + A( L, I )*B( L, J )
258   100             CONTINUE
259                   IF( BETA.EQ.ZERO )THEN
260                      C( I, J ) = ALPHA*TEMP
261                   ELSE
262                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
263                   END IF
264   110          CONTINUE
265   120       CONTINUE
266          END IF
267       ELSE
268          IF( NOTA )THEN
269 *
270 *           Form  C := alpha*A*B' + beta*C
271 *
272             DO 170, J = 1, N
273                IF( BETA.EQ.ZERO )THEN
274                   DO 130, I = 1, M
275                      C( I, J ) = ZERO
276   130             CONTINUE
277                ELSE IF( BETA.NE.ONE )THEN
278                   DO 140, I = 1, M
279                      C( I, J ) = BETA*C( I, J )
280   140             CONTINUE
281                END IF
282                DO 160, L = 1, K
283                   IF( B( J, L ).NE.ZERO )THEN
284                      TEMP = ALPHA*B( J, L )
285                      DO 150, I = 1, M
286                         C( I, J ) = C( I, J ) + TEMP*A( I, L )
287   150                CONTINUE
288                   END IF
289   160          CONTINUE
290   170       CONTINUE
291          ELSE
292 *
293 *           Form  C := alpha*A'*B' + beta*C
294 *
295             DO 200, J = 1, N
296                DO 190, I = 1, M
297                   TEMP = ZERO
298                   DO 180, L = 1, K
299                      TEMP = TEMP + A( L, I )*B( J, L )
300   180             CONTINUE
301                   IF( BETA.EQ.ZERO )THEN
302                      C( I, J ) = ALPHA*TEMP
303                   ELSE
304                      C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
305                   END IF
306   190          CONTINUE
307   200       CONTINUE
308          END IF
309       END IF
310 *
311       RETURN
312 *
313 *     End of DGEMM .
314 *
315       END
316
317       SUBROUTINE OMXUNSAFEDGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
318      $                   BETA, Y, INCY )
319 *     .. Scalar Arguments ..
320       DOUBLE PRECISION   ALPHA, BETA
321       INTEGER            INCX, INCY, LDA, M, N
322       CHARACTER*1        TRANS
323 *     .. Array Arguments ..
324       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
325 *     ..
326 *
327 *  Purpose
328 *  =======
329 *
330 *  DGEMV  performs one of the matrix-vector operations
331 *
332 *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
333 *
334 *  where alpha and beta are scalars, x and y are vectors and A is an
335 *  m by n matrix.
336 *
337 *  Parameters
338 *  ==========
339 *
340 *  TRANS  - CHARACTER*1.
341 *           On entry, TRANS specifies the operation to be performed as
342 *           follows:
343 *
344 *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
345 *
346 *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
347 *
348 *              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
349 *
350 *           Unchanged on exit.
351 *
352 *  M      - INTEGER.
353 *           On entry, M specifies the number of rows of the matrix A.
354 *           M must be at least zero.
355 *           Unchanged on exit.
356 *
357 *  N      - INTEGER.
358 *           On entry, N specifies the number of columns of the matrix A.
359 *           N must be at least zero.
360 *           Unchanged on exit.
361 *
362 *  ALPHA  - DOUBLE PRECISION.
363 *           On entry, ALPHA specifies the scalar alpha.
364 *           Unchanged on exit.
365 *
366 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
367 *           Before entry, the leading m by n part of the array A must
368 *           contain the matrix of coefficients.
369 *           Unchanged on exit.
370 *
371 *  LDA    - INTEGER.
372 *           On entry, LDA specifies the first dimension of A as declared
373 *           in the calling (sub) program. LDA must be at least
374 *           max( 1, m ).
375 *           Unchanged on exit.
376 *
377 *  X      - DOUBLE PRECISION array of DIMENSION at least
378 *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
379 *           and at least
380 *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
381 *           Before entry, the incremented array X must contain the
382 *           vector x.
383 *           Unchanged on exit.
384 *
385 *  INCX   - INTEGER.
386 *           On entry, INCX specifies the increment for the elements of
387 *           X. INCX must not be zero.
388 *           Unchanged on exit.
389 *
390 *  BETA   - DOUBLE PRECISION.
391 *           On entry, BETA specifies the scalar beta. When BETA is
392 *           supplied as zero then Y need not be set on input.
393 *           Unchanged on exit.
394 *
395 *  Y      - DOUBLE PRECISION array of DIMENSION at least
396 *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
397 *           and at least
398 *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
399 *           Before entry with BETA non-zero, the incremented array Y
400 *           must contain the vector y. On exit, Y is overwritten by the
401 *           updated vector y.
402 *
403 *  INCY   - INTEGER.
404 *           On entry, INCY specifies the increment for the elements of
405 *           Y. INCY must not be zero.
406 *           Unchanged on exit.
407 *
408 *
409 *  Level 2 Blas routine.
410 *
411 *  -- Written on 22-October-1986.
412 *     Jack Dongarra, Argonne National Lab.
413 *     Jeremy Du Croz, Nag Central Office.
414 *     Sven Hammarling, Nag Central Office.
415 *     Richard Hanson, Sandia National Labs.
416 *
417 *
418 *     .. Parameters ..
419       DOUBLE PRECISION   ONE         , ZERO
420       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
421 *     .. Local Scalars ..
422       DOUBLE PRECISION   TEMP
423       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
424 *     .. External Functions ..
425       LOGICAL            LSAME
426       EXTERNAL           LSAME
427 *     .. External Subroutines ..
428       EXTERNAL           XERBLA
429 *     .. Intrinsic Functions ..
430       INTRINSIC          MAX
431 *     ..
432 *     .. Executable Statements ..
433 *
434 *     Test the input parameters.
435 *
436       INFO = 0
437       IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
438      $         .NOT.LSAME( TRANS, 'T' ).AND.
439      $         .NOT.LSAME( TRANS, 'C' )      )THEN
440          INFO = 1
441       ELSE IF( M.LT.0 )THEN
442          INFO = 2
443       ELSE IF( N.LT.0 )THEN
444          INFO = 3
445       ELSE IF( LDA.LT.MAX( 1, M ) )THEN
446          INFO = 6
447       ELSE IF( INCX.EQ.0 )THEN
448          INFO = 8
449       ELSE IF( INCY.EQ.0 )THEN
450          INFO = 11
451       END IF
452       IF( INFO.NE.0 )THEN
453          CALL XERBLA( 'OMXUNSAFEDGEMV ', INFO )
454          RETURN
455       END IF
456 *
457 *     Quick return if possible.
458 *
459       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
460      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
461      $   RETURN
462 *
463 *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
464 *     up the start points in  X  and  Y.
465 *
466       IF( LSAME( TRANS, 'N' ) )THEN
467          LENX = N
468          LENY = M
469       ELSE
470          LENX = M
471          LENY = N
472       END IF
473       IF( INCX.GT.0 )THEN
474          KX = 1
475       ELSE
476          KX = 1 - ( LENX - 1 )*INCX
477       END IF
478       IF( INCY.GT.0 )THEN
479          KY = 1
480       ELSE
481          KY = 1 - ( LENY - 1 )*INCY
482       END IF
483 *
484 *     Start the operations. In this version the elements of A are
485 *     accessed sequentially with one pass through A.
486 *
487 *     First form  y := beta*y.
488 *
489       IF( BETA.NE.ONE )THEN
490          IF( INCY.EQ.1 )THEN
491             IF( BETA.EQ.ZERO )THEN
492                DO 10, I = 1, LENY
493                   Y( I ) = ZERO
494    10          CONTINUE
495             ELSE
496                DO 20, I = 1, LENY
497                   Y( I ) = BETA*Y( I )
498    20          CONTINUE
499             END IF
500          ELSE
501             IY = KY
502             IF( BETA.EQ.ZERO )THEN
503                DO 30, I = 1, LENY
504                   Y( IY ) = ZERO
505                   IY      = IY   + INCY
506    30          CONTINUE
507             ELSE
508                DO 40, I = 1, LENY
509                   Y( IY ) = BETA*Y( IY )
510                   IY      = IY           + INCY
511    40          CONTINUE
512             END IF
513          END IF
514       END IF
515       IF( ALPHA.EQ.ZERO )
516      $   RETURN
517       IF( LSAME( TRANS, 'N' ) )THEN
518 *
519 *        Form  y := alpha*A*x + y.
520 *
521          JX = KX
522         IF( INCY.EQ.1 )THEN
523             DO 60, J = 1, N
524                IF( X( JX ).NE.ZERO )THEN
525                   TEMP = ALPHA*X( JX )
526                   DO 50, I = 1, M
527                      Y( I ) = Y( I ) + TEMP*A( I, J )
528    50             CONTINUE
529                END IF
530                JX = JX + INCX
531    60       CONTINUE
532          ELSE
533             DO 80, J = 1, N
534                IF( X( JX ).NE.ZERO )THEN
535                   TEMP = ALPHA*X( JX )
536                   IY   = KY
537                   DO 70, I = 1, M
538                      Y( IY ) = Y( IY ) + TEMP*A( I, J )
539                      IY      = IY      + INCY
540    70             CONTINUE
541                END IF
542                JX = JX + INCX
543    80       CONTINUE
544          END IF
545       ELSE
546 *
547 *        Form  y := alpha*A'*x + y.
548 *
549          JY = KY
550          IF( INCX.EQ.1 )THEN
551             DO 100, J = 1, N
552                TEMP = ZERO
553                DO 90, I = 1, M
554                   TEMP = TEMP + A( I, J )*X( I )
555    90          CONTINUE
556                Y( JY ) = Y( JY ) + ALPHA*TEMP
557                JY      = JY      + INCY
558   100       CONTINUE
559          ELSE
560             DO 120, J = 1, N
561                TEMP = ZERO
562                IX   = KX
563         DO 110, I = 1, M
564                   TEMP = TEMP + A( I, J )*X( IX )
565                   IX   = IX   + INCX
566   110          CONTINUE
567                Y( JY ) = Y( JY ) + ALPHA*TEMP
568                JY      = JY      + INCY
569   120       CONTINUE
570          END IF
571       END IF
572 *
573       RETURN
574 *
575 *     End of DGEMV .
576 *
577       END
578