Added some missing headers reported by /usr/lib/build/checks-data/check_gcc_output
[scilab.git] / scilab / modules / scicos_blocks / src / c / dmmul.c
1 /*  Scicos
2 *
3 *  Copyright (C) INRIA - METALAU Project <scicos@inria.fr>
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18 *
19 * See the file ./license.txt
20 */
21 /*--------------------------------------------------------------------------*/ 
22 #include <stdio.h> /* printf */
23 #include <string.h>
24 #ifdef UNIX
25   #include <unistd.h>
26   #include <sys/socket.h>
27 #endif
28 #include "scicos_math.h"
29 #include "dynlib_scicos_blocks.h"
30 /*--------------------------------------------------------------------------*/ 
31 /* Table of constant values */
32 static double c_b4 = 1.;
33 static double c_b5 = 0.;
34 /*--------------------------------------------------------------------------*/ 
35 static int dgemm(char *transa, char *transb, int *m, int *n, int *k, double *alpha, 
36                   double *a, int *lda, double *b, int *ldb,double *beta, double *c, 
37                   int *ldc);
38 static long int lsame(char *ca, char *cb);
39 static int xerbla(char *srname, int *info);
40 /*--------------------------------------------------------------------------*/ 
41 SCICOS_BLOCKS_IMPEXP int dmmul(double *a, int *na, double *b, int *nb, double *c__, 
42           int *nc, int *l, int *m, int *n)
43 {
44   int a_dim1 = 0, a_offset = 0, b_dim1 = 0, b_offset = 0, c_dim1 = 0, c_offset = 0;
45     
46 /*     PURPOSE */
47 /*        computes the matrix product C = A * B */
48 /*            C   =   A   *   B */
49 /*          (l,n)   (l,m) * (m,n) */
50
51 /*     PARAMETERS */
52 /*        input */
53 /*        ----- */
54 /*        A : (double) array (l, m) with leading dim na */
55
56 /*        B : (double) array (m, n) with leading dim nb */
57
58 /*        na, nb, nc, l, m, n : integers */
59
60 /*        output */
61 /*        ------ */
62 /*        C : (double) array (l, n) with leading dim nc */
63
64 /*     NOTE */
65 /*        (original version substituted by a call to the blas dgemm) */
66     /* Parameter adjustments */
67     a_dim1 = *na;
68     a_offset = 1 + a_dim1 * 1;
69     a -= a_offset;
70     c_dim1 = *nc;
71     c_offset = 1 + c_dim1 * 1;
72     c__ -= c_offset;
73     b_dim1 = *nb;
74     b_offset = 1 + b_dim1 * 1;
75     b -= b_offset;
76
77     /* Function Body */
78     dgemm("n", "n", l, n, m, &c_b4, &a[a_offset], na, &b[b_offset], nb, &
79             c_b5, &c__[c_offset], nc);
80         return 0;
81 } /* dmmul */
82 /*--------------------------------------------------------------------------*/ 
83 static int dgemm(char *transa, char *transb, int *m, int *n, int *k, double *alpha, 
84           double *a, int *lda, double *b, int *ldb,double *beta, double *c, 
85           int *ldc)
86 {
87   /* System generated locals */
88   int i__1, i__2, i__3;
89   /* Local variables */
90   static int info = 0;
91   static long int nota = 0, notb = 0;
92   static double temp = 0.;
93   static int i = 0, j = 0, l = 0, ncola = 0;
94   static int nrowa = 0, nrowb = 0;  
95
96   
97   /*     .. Scalar Arguments .. */
98   /*     .. Array Arguments .. */
99   /*     .. */
100   
101   /*  Purpose */
102   /*  ======= */
103   
104   /*  DGEMM  performs one of the matrix-matrix operations */
105   
106   /*     C := alpha*op( A )*op( B ) + beta*C, */
107   
108   /*  where  op( X ) is one of */
109   
110   /*     op( X ) = X   or   op( X ) = X', */
111   
112   /*  alpha and beta are scalars, and A, B and C are matrices, with op( A ) */
113   /*  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix. */
114   
115   /*  Parameters */
116   /*  ========== */
117   
118   /*  TRANSA - CHARACTER*1. */
119   /*           On entry, TRANSA specifies the form of op( A ) to be used in */
120   /*           the matrix multiplication as follows: */
121   
122   /*              TRANSA = 'N' or 'n',  op( A ) = A. */
123   
124   /*              TRANSA = 'T' or 't',  op( A ) = A'. */
125   
126   /*              TRANSA = 'C' or 'c',  op( A ) = A'. */
127   
128   /*           Unchanged on exit. */
129   
130   /*  TRANSB - CHARACTER*1. */
131   /*           On entry, TRANSB specifies the form of op( B ) to be used in */
132   /*           the matrix multiplication as follows: */
133   
134   /*              TRANSB = 'N' or 'n',  op( B ) = B. */
135   
136   /*              TRANSB = 'T' or 't',  op( B ) = B'. */
137   
138   /*              TRANSB = 'C' or 'c',  op( B ) = B'. */
139   
140   /*           Unchanged on exit. */
141   
142   /*  M      - INTEGER. */
143   /*           On entry,  M  specifies  the number  of rows  of the  matrix */
144   /*           op( A )  and of the  matrix  C.  M  must  be at least  zero. */
145   /*           Unchanged on exit. */
146   
147   /*  N      - INTEGER. */
148   /*           On entry,  N  specifies the number  of columns of the matrix */
149   /*           op( B ) and the number of columns of the matrix C. N must be */
150   /*           at least zero. */
151   /*           Unchanged on exit. */
152   
153   /*  K      - INTEGER. */
154   /*           On entry,  K  specifies  the number of columns of the matrix */
155   /*           op( A ) and the number of rows of the matrix op( B ). K must */
156   /*           be at least  zero. */
157   /*           Unchanged on exit. */
158   
159   /*  ALPHA  - DOUBLE PRECISION. */
160   /*           On entry, ALPHA specifies the scalar alpha. */
161   /*           Unchanged on exit. */
162   
163   /*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is */
164   /*           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise. */
165   /*           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k */
166   /*           part of the array  A  must contain the matrix  A,  otherwise */
167   /*           the leading  k by m  part of the array  A  must contain  the */
168   /*           matrix A. */
169   /*           Unchanged on exit. */
170   
171   /*  LDA    - INTEGER. */
172   /*           On entry, LDA specifies the first dimension of A as declared */
173   /*           in the calling (sub) program. When  TRANSA = 'N' or 'n' then */
174   /*           LDA must be at least  max( 1, m ), otherwise  LDA must be at */
175   /*           least  max( 1, k ). */
176   /*           Unchanged on exit. */
177   
178   /*  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is */
179   /*           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise. */
180   /*           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n */
181   /*           part of the array  B  must contain the matrix  B,  otherwise */
182   /*           the leading  n by k  part of the array  B  must contain  the */
183   /*           matrix B. */
184   /*           Unchanged on exit. */
185   
186   /*  LDB    - INTEGER. */
187   /*           On entry, LDB specifies the first dimension of B as declared */
188   /*           in the calling (sub) program. When  TRANSB = 'N' or 'n' then */
189   /*           LDB must be at least  max( 1, k ), otherwise  LDB must be at */
190   /*           least  max( 1, n ). */
191   /*           Unchanged on exit. */
192   
193   /*  BETA   - DOUBLE PRECISION. */
194   /*           On entry,  BETA  specifies the scalar  beta.  When  BETA  is */
195   /*           supplied as zero then C need not be set on input. */
196   /*           Unchanged on exit. */
197   
198   /*  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ). */
199   /*           Before entry, the leading  m by n  part of the array  C must */
200   /*           contain the matrix  C,  except when  beta  is zero, in which */
201   /*           case C need not be set on entry. */
202   /*           On exit, the array  C  is overwritten by the  m by n  matrix */
203   /*           ( alpha*op( A )*op( B ) + beta*C ). */
204   
205   /*  LDC    - INTEGER. */
206   /*           On entry, LDC specifies the first dimension of C as declared */
207   /*           in  the  calling  (sub)  program.   LDC  must  be  at  least */
208   /*           max( 1, m ). */
209   /*           Unchanged on exit. */
210   
211   
212   /*  Level 3 Blas routine. */
213   
214   /*  -- Written on 8-February-1989. */
215   /*     Jack Dongarra, Argonne National Laboratory. */
216   /*     Iain Duff, AERE Harwell. */
217   /*     Jeremy Du Croz, Numerical Algorithms Group Ltd. */
218   /*     Sven Hammarling, Numerical Algorithms Group Ltd. */
219   
220   
221   /*     .. External Functions .. */
222   /*     .. External Subroutines .. */
223   /*     .. Intrinsic Functions .. */
224   /*     .. Local Scalars .. */
225   /*     .. Parameters .. */
226   /*     .. */
227   /*     .. Executable Statements .. */
228   
229   /*     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not */
230   /*     transposed and set  NROWA, NCOLA and  NROWB  as the number of rows */
231   /*     and  columns of  A  and the  number of  rows  of  B  respectively. */
232   
233 #define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]
234 #define B(I,J) b[(I)-1 + ((J)-1)* ( *ldb)]
235 #define C(I,J) c[(I)-1 + ((J)-1)* ( *ldc)]
236   
237   nota = lsame(transa, "N");
238   notb = lsame(transb, "N");
239   if (nota) {
240     nrowa = *m;
241     ncola = *k;
242   } else {
243     nrowa = *k;
244     ncola = *m;
245   }
246   if (notb) {
247     nrowb = *k;
248   } else {
249     nrowb = *n;
250   }
251   
252   /*     Test the input parameters. */
253   
254   info = 0;
255   if (! nota && ! lsame(transa, "C") && ! lsame(transa, "T")) {
256     info = 1;
257   } else if (! notb && ! lsame(transb, "C") && ! lsame(transb,"T")) {
258     info = 2;
259   } else if (*m < 0) {
260     info = 3;
261   } else if (*n < 0) {
262     info = 4;
263   } else if (*k < 0) {
264     info = 5;
265   } else if (*lda < max(1,nrowa)) {
266     info = 8;
267   } else if (*ldb < max(1,nrowb)) {
268     info = 10;
269   } else if (*ldc < max(1,*m)) {
270     info = 13;
271   }
272   if (info != 0) {
273     xerbla("DGEMM ", &info);
274     return 0;
275   }
276   
277   /*     Quick return if possible. */
278   
279   if (*m == 0 || *n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.) {
280     return 0;
281   }
282   
283   /*     And if  alpha.eq.zero. */
284   
285   if (*alpha == 0.) {
286     if (*beta == 0.) {
287       i__1 = *n;
288       for (j = 1; j <= *n; ++j) {
289         i__2 = *m;
290         for (i = 1; i <= *m; ++i) {
291           C(i,j) = 0.;
292           /* L10: */
293         }
294         /* L20: */
295       }
296     } else {
297       i__1 = *n;
298       for (j = 1; j <= *n; ++j) {
299         i__2 = *m;
300         for (i = 1; i <= *m; ++i) {
301           C(i,j) = *beta * C(i,j);
302           /* L30: */
303         }
304         /* L40: */
305       }
306     }
307     return 0;
308   }
309   
310   /*     Start the operations. */
311   
312   if (notb) {
313     if (nota) {
314       
315       /*           Form  C := alpha*A*B + beta*C. */
316       
317       i__1 = *n;
318       for (j = 1; j <= *n; ++j) {
319         if (*beta == 0.) {
320           i__2 = *m;
321           for (i = 1; i <= *m; ++i) {
322             C(i,j) = 0.;
323             /* L50: */
324           }
325         } else if (*beta != 1.) {
326           i__2 = *m;
327           for (i = 1; i <= *m; ++i) {
328             C(i,j) = *beta * C(i,j);
329             /* L60: */
330           }
331         }
332         i__2 = *k;
333         for (l = 1; l <= *k; ++l) {
334           if (B(l,j) != 0.) {
335             temp = *alpha * B(l,j);
336             i__3 = *m;
337             for (i = 1; i <= *m; ++i) {
338               C(i,j) += temp * A(i,l);
339               /* L70: */
340             }
341           }
342           /* L80: */
343         }
344         /* L90: */
345       }
346     } else {
347       
348       /*           Form  C := alpha*A'*B + beta*C */
349       
350       i__1 = *n;
351       for (j = 1; j <= *n; ++j) {
352         i__2 = *m;
353         for (i = 1; i <= *m; ++i) {
354           temp = 0.;
355           i__3 = *k;
356           for (l = 1; l <= *k; ++l) {
357             temp += A(l,i) * B(l,j);
358             /* L100: */
359           }
360           if (*beta == 0.) {
361             C(i,j) = *alpha * temp;
362           } else {
363             C(i,j) = *alpha * temp + *beta * C(i,j);
364           }
365           /* L110: */
366         }
367         /* L120: */
368       }
369     }
370   } else {
371     if (nota) {
372       
373       /*           Form  C := alpha*A*B' + beta*C */
374       
375       i__1 = *n;
376       for (j = 1; j <= *n; ++j) {
377         if (*beta == 0.) {
378           i__2 = *m;
379           for (i = 1; i <= *m; ++i) {
380             C(i,j) = 0.;
381             /* L130: */
382           }
383         } else if (*beta != 1.) {
384           i__2 = *m;
385           for (i = 1; i <= *m; ++i) {
386             C(i,j) = *beta * C(i,j);
387             /* L140: */
388           }
389         }
390         i__2 = *k;
391         for (l = 1; l <= *k; ++l) {
392           if (B(j,l) != 0.) {
393             temp = *alpha * B(j,l);
394             i__3 = *m;
395             for (i = 1; i <= *m; ++i) {
396               C(i,j) += temp * A(i,l);
397               /* L150: */
398             }
399           }
400           /* L160: */
401         }
402         /* L170: */
403       }
404     } else {
405       
406       /*           Form  C := alpha*A'*B' + beta*C */
407       
408       i__1 = *n;
409       for (j = 1; j <= *n; ++j) {
410         i__2 = *m;
411         for (i = 1; i <= *m; ++i) {
412           temp = 0.;
413           i__3 = *k;
414           for (l = 1; l <= *k; ++l) {
415             temp += A(l,i) * B(j,l);
416             /* L180: */
417           }
418           if (*beta == 0.) {
419             C(i,j) = *alpha * temp;
420           } else {
421             C(i,j) = *alpha * temp + *beta * C(i,j);
422           }
423           /* L190: */
424         }
425         /* L200: */
426       }
427     }
428   }
429   
430   return 0;
431   
432   /*     End of DGEMM . */
433   
434 } /* dgemm */
435 /*--------------------------------------------------------------------------*/ 
436 static long int lsame(char *ca, char *cb)
437 {
438   /* System generated locals */
439   long int ret_val = 0;
440
441   /* Local variables */
442   static int inta = 0, intb = 0, zcode = 0;
443   
444
445   /*  -- LAPACK auxiliary routine (version 2.0) -- */
446   /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
447   /*     Courant Institute, Argonne National Lab, and Rice University */
448   /*     January 31, 1994 */
449   
450   /*     .. Scalar Arguments .. */
451   /*     .. */
452   
453   /*  Purpose */
454   /*  ======= */
455   
456   /*  LSAME returns .TRUE. if CA is the same letter as CB regardless of */
457   /*  case. */
458   
459   /*  Arguments */
460   /*  ========= */
461   
462   /*  CA      (input) CHARACTER*1 */
463   /*  CB      (input) CHARACTER*1 */
464   /*          CA and CB specify the single characters to be compared. */
465   
466   /* ===================================================================== */
467   
468   /*     .. Intrinsic Functions .. */
469   /*     .. */
470   /*     .. Local Scalars .. */
471   /*     .. */
472   /*     .. Executable Statements .. */
473   
474   /*     Test if the characters are equal */
475   
476   ret_val = *(unsigned char *)ca == *(unsigned char *)cb;
477   if (ret_val) {
478     return ret_val;
479   }
480   
481   /*     Now test for equivalence if both characters are alphabetic. */
482   
483   zcode = 'Z';
484   
485   /*     Use 'Z' rather than 'A' so that ASCII can be detected on Prime */
486   /*     machines, on which ICHAR returns a value with bit 8 set. */
487   /*     ICHAR('A') on Prime machines returns 193 which is the same as */
488   /*     ICHAR('A') on an EBCDIC machine. */
489   
490   inta = *(unsigned char *)ca;
491   intb = *(unsigned char *)cb;
492   
493   if (zcode == 90 || zcode == 122) {
494     
495     /*        ASCII is assumed - ZCODE is the ASCII code of either lower or */
496     /*        upper case 'Z'. */
497     
498     if (inta >= 97 && inta <= 122) {
499       inta += -32;
500     }
501     if (intb >= 97 && intb <= 122) {
502       intb += -32;
503     }
504     
505   } else if (zcode == 233 || zcode == 169) {
506     
507     /*        EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or */
508     /*        upper case 'Z'. */
509     
510     if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta 
511         >= 162 && inta <= 169) {
512       inta += 64;
513     }
514     if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb 
515         >= 162 && intb <= 169) {
516       intb += 64;
517     }
518     
519   } else if (zcode == 218 || zcode == 250) {
520     
521     /*        ASCII is assumed, on Prime machines - ZCODE is the ASCII code */
522     /*        plus 128 of either lower or upper case 'Z'. */
523     
524     if (inta >= 225 && inta <= 250) {
525       inta += -32;
526     }
527     if (intb >= 225 && intb <= 250) {
528       intb += -32;
529     }
530   }
531   ret_val = inta == intb;
532   
533   /*     RETURN */
534   
535   /*     End of LSAME */
536   
537   return ret_val;
538 } /* lsame */
539 /*--------------------------------------------------------------------------*/ 
540 static int xerbla(char *srname, int *info)
541 {
542   
543   /*  -- LAPACK auxiliary routine (version 3.0) -- */
544   /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
545   /*     Courant Institute, Argonne National Lab, and Rice University */
546   /*     September 30, 1994 */
547   
548   /*     .. Scalar Arguments .. */
549   /*     .. */
550   
551   /*  Purpose */
552   /*  ======= */
553   
554   /*  XERBLA  is an error handler for the LAPACK routines. */
555   /*  It is called by an LAPACK routine if an input parameter has an */
556   /*  invalid value.  A message is printed and execution stops. */
557   
558   /*  Installers may consider modifying the STOP statement in order to */
559   /*  call system-specific exception-handling facilities. */
560   
561   /*  Arguments */
562   /*  ========= */
563   
564   /*  SRNAME  (input) CHARACTER*6 */
565   /*          The name of the routine which called XERBLA. */
566   
567   /*  INFO    (input) INTEGER */
568   /*          The position of the invalid parameter in the parameter list */
569   /*          of the calling routine. */
570   
571   /* ===================================================================== */
572   
573   printf("** On entry to %6s, parameter number %2i had an illegal value\n",
574          srname, *info);
575   
576   
577   /*     End of XERBLA */
578   return 0; 
579 } /* xerbla */
580 /*--------------------------------------------------------------------------*/