~ubuntu-branches/ubuntu/trusty/igraph/trusty-proposed

« back to all changes in this revision

Viewing changes to src/lapack/dtrmm.c

  • Committer: Package Import Robot
  • Author(s): Mathieu Malaterre
  • Date: 2013-05-27 14:01:54 UTC
  • mfrom: (4.1.1 experimental)
  • Revision ID: package-import@ubuntu.com-20130527140154-oxwwmr0gj3kdy4ol
Tags: 0.6.5-2
Upload to sid

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/*  -- translated by f2c (version 20100827).
 
2
   You must link the resulting object file with libf2c:
 
3
        on Microsoft Windows system, link with libf2c.lib;
 
4
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
 
5
        or, if you install libf2c.a in a standard place, with -lf2c -lm
 
6
        -- in that order, at the end of the command line, as in
 
7
                cc *.o -lf2c -lm
 
8
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
 
9
 
 
10
                http://www.netlib.org/f2c/libf2c.zip
 
11
*/
 
12
 
 
13
#include "f2c.h"
 
14
 
 
15
/* Subroutine */ int igraphdtrmm_(char *side, char *uplo, char *transa, char *diag, 
 
16
        integer *m, integer *n, doublereal *alpha, doublereal *a, integer *
 
17
        lda, doublereal *b, integer *ldb)
 
18
{
 
19
    /* System generated locals */
 
20
    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
 
21
 
 
22
    /* Local variables */
 
23
    integer i__, j, k, info;
 
24
    doublereal temp;
 
25
    logical lside;
 
26
    extern logical igraphlsame_(char *, char *);
 
27
    integer nrowa;
 
28
    logical upper;
 
29
    extern /* Subroutine */ int igraphxerbla_(char *, integer *, ftnlen);
 
30
    logical nounit;
 
31
 
 
32
 
 
33
/*  Purpose   
 
34
    =======   
 
35
 
 
36
    DTRMM  performs one of the matrix-matrix operations   
 
37
 
 
38
       B := alpha*op( A )*B,   or   B := alpha*B*op( A ),   
 
39
 
 
40
    where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or   
 
41
    non-unit,  upper or lower triangular matrix  and  op( A )  is one  of   
 
42
 
 
43
       op( A ) = A   or   op( A ) = A**T.   
 
44
 
 
45
    Arguments   
 
46
    ==========   
 
47
 
 
48
    SIDE   - CHARACTER*1.   
 
49
             On entry,  SIDE specifies whether  op( A ) multiplies B from   
 
50
             the left or right as follows:   
 
51
 
 
52
                SIDE = 'L' or 'l'   B := alpha*op( A )*B.   
 
53
 
 
54
                SIDE = 'R' or 'r'   B := alpha*B*op( A ).   
 
55
 
 
56
             Unchanged on exit.   
 
57
 
 
58
    UPLO   - CHARACTER*1.   
 
59
             On entry, UPLO specifies whether the matrix A is an upper or   
 
60
             lower triangular matrix as follows:   
 
61
 
 
62
                UPLO = 'U' or 'u'   A is an upper triangular matrix.   
 
63
 
 
64
                UPLO = 'L' or 'l'   A is a lower triangular matrix.   
 
65
 
 
66
             Unchanged on exit.   
 
67
 
 
68
    TRANSA - CHARACTER*1.   
 
69
             On entry, TRANSA specifies the form of op( A ) to be used in   
 
70
             the matrix multiplication as follows:   
 
71
 
 
72
                TRANSA = 'N' or 'n'   op( A ) = A.   
 
73
 
 
74
                TRANSA = 'T' or 't'   op( A ) = A**T.   
 
75
 
 
76
                TRANSA = 'C' or 'c'   op( A ) = A**T.   
 
77
 
 
78
             Unchanged on exit.   
 
79
 
 
80
    DIAG   - CHARACTER*1.   
 
81
             On entry, DIAG specifies whether or not A is unit triangular   
 
82
             as follows:   
 
83
 
 
84
                DIAG = 'U' or 'u'   A is assumed to be unit triangular.   
 
85
 
 
86
                DIAG = 'N' or 'n'   A is not assumed to be unit   
 
87
                                    triangular.   
 
88
 
 
89
             Unchanged on exit.   
 
90
 
 
91
    M      - INTEGER.   
 
92
             On entry, M specifies the number of rows of B. M must be at   
 
93
             least zero.   
 
94
             Unchanged on exit.   
 
95
 
 
96
    N      - INTEGER.   
 
97
             On entry, N specifies the number of columns of B.  N must be   
 
98
             at least zero.   
 
99
             Unchanged on exit.   
 
100
 
 
101
    ALPHA  - DOUBLE PRECISION.   
 
102
             On entry,  ALPHA specifies the scalar  alpha. When  alpha is   
 
103
             zero then  A is not referenced and  B need not be set before   
 
104
             entry.   
 
105
             Unchanged on exit.   
 
106
 
 
107
    A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m   
 
108
             when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.   
 
109
             Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k   
 
110
             upper triangular part of the array  A must contain the upper   
 
111
             triangular matrix  and the strictly lower triangular part of   
 
112
             A is not referenced.   
 
113
             Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k   
 
114
             lower triangular part of the array  A must contain the lower   
 
115
             triangular matrix  and the strictly upper triangular part of   
 
116
             A is not referenced.   
 
117
             Note that when  DIAG = 'U' or 'u',  the diagonal elements of   
 
118
             A  are not referenced either,  but are assumed to be  unity.   
 
119
             Unchanged on exit.   
 
120
 
 
121
    LDA    - INTEGER.   
 
122
             On entry, LDA specifies the first dimension of A as declared   
 
123
             in the calling (sub) program.  When  SIDE = 'L' or 'l'  then   
 
124
             LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'   
 
125
             then LDA must be at least max( 1, n ).   
 
126
             Unchanged on exit.   
 
127
 
 
128
    B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).   
 
129
             Before entry,  the leading  m by n part of the array  B must   
 
130
             contain the matrix  B,  and  on exit  is overwritten  by the   
 
131
             transformed matrix.   
 
132
 
 
133
    LDB    - INTEGER.   
 
134
             On entry, LDB specifies the first dimension of B as declared   
 
135
             in  the  calling  (sub)  program.   LDB  must  be  at  least   
 
136
             max( 1, m ).   
 
137
             Unchanged on exit.   
 
138
 
 
139
    Further Details   
 
140
    ===============   
 
141
 
 
142
    Level 3 Blas routine.   
 
143
 
 
144
    -- Written on 8-February-1989.   
 
145
       Jack Dongarra, Argonne National Laboratory.   
 
146
       Iain Duff, AERE Harwell.   
 
147
       Jeremy Du Croz, Numerical Algorithms Group Ltd.   
 
148
       Sven Hammarling, Numerical Algorithms Group Ltd.   
 
149
 
 
150
    =====================================================================   
 
151
 
 
152
 
 
153
       Test the input parameters.   
 
154
 
 
155
       Parameter adjustments */
 
156
    a_dim1 = *lda;
 
157
    a_offset = 1 + a_dim1;
 
158
    a -= a_offset;
 
159
    b_dim1 = *ldb;
 
160
    b_offset = 1 + b_dim1;
 
161
    b -= b_offset;
 
162
 
 
163
    /* Function Body */
 
164
    lside = igraphlsame_(side, "L");
 
165
    if (lside) {
 
166
        nrowa = *m;
 
167
    } else {
 
168
        nrowa = *n;
 
169
    }
 
170
    nounit = igraphlsame_(diag, "N");
 
171
    upper = igraphlsame_(uplo, "U");
 
172
 
 
173
    info = 0;
 
174
    if (! lside && ! igraphlsame_(side, "R")) {
 
175
        info = 1;
 
176
    } else if (! upper && ! igraphlsame_(uplo, "L")) {
 
177
        info = 2;
 
178
    } else if (! igraphlsame_(transa, "N") && ! igraphlsame_(transa,
 
179
             "T") && ! igraphlsame_(transa, "C")) {
 
180
        info = 3;
 
181
    } else if (! igraphlsame_(diag, "U") && ! igraphlsame_(diag, 
 
182
            "N")) {
 
183
        info = 4;
 
184
    } else if (*m < 0) {
 
185
        info = 5;
 
186
    } else if (*n < 0) {
 
187
        info = 6;
 
188
    } else if (*lda < max(1,nrowa)) {
 
189
        info = 9;
 
190
    } else if (*ldb < max(1,*m)) {
 
191
        info = 11;
 
192
    }
 
193
    if (info != 0) {
 
194
        igraphxerbla_("DTRMM ", &info, (ftnlen)6);
 
195
        return 0;
 
196
    }
 
197
 
 
198
/*     Quick return if possible. */
 
199
 
 
200
    if (*m == 0 || *n == 0) {
 
201
        return 0;
 
202
    }
 
203
 
 
204
/*     And when  alpha.eq.zero. */
 
205
 
 
206
    if (*alpha == 0.) {
 
207
        i__1 = *n;
 
208
        for (j = 1; j <= i__1; ++j) {
 
209
            i__2 = *m;
 
210
            for (i__ = 1; i__ <= i__2; ++i__) {
 
211
                b[i__ + j * b_dim1] = 0.;
 
212
/* L10: */
 
213
            }
 
214
/* L20: */
 
215
        }
 
216
        return 0;
 
217
    }
 
218
 
 
219
/*     Start the operations. */
 
220
 
 
221
    if (lside) {
 
222
        if (igraphlsame_(transa, "N")) {
 
223
 
 
224
/*           Form  B := alpha*A*B. */
 
225
 
 
226
            if (upper) {
 
227
                i__1 = *n;
 
228
                for (j = 1; j <= i__1; ++j) {
 
229
                    i__2 = *m;
 
230
                    for (k = 1; k <= i__2; ++k) {
 
231
                        if (b[k + j * b_dim1] != 0.) {
 
232
                            temp = *alpha * b[k + j * b_dim1];
 
233
                            i__3 = k - 1;
 
234
                            for (i__ = 1; i__ <= i__3; ++i__) {
 
235
                                b[i__ + j * b_dim1] += temp * a[i__ + k * 
 
236
                                        a_dim1];
 
237
/* L30: */
 
238
                            }
 
239
                            if (nounit) {
 
240
                                temp *= a[k + k * a_dim1];
 
241
                            }
 
242
                            b[k + j * b_dim1] = temp;
 
243
                        }
 
244
/* L40: */
 
245
                    }
 
246
/* L50: */
 
247
                }
 
248
            } else {
 
249
                i__1 = *n;
 
250
                for (j = 1; j <= i__1; ++j) {
 
251
                    for (k = *m; k >= 1; --k) {
 
252
                        if (b[k + j * b_dim1] != 0.) {
 
253
                            temp = *alpha * b[k + j * b_dim1];
 
254
                            b[k + j * b_dim1] = temp;
 
255
                            if (nounit) {
 
256
                                b[k + j * b_dim1] *= a[k + k * a_dim1];
 
257
                            }
 
258
                            i__2 = *m;
 
259
                            for (i__ = k + 1; i__ <= i__2; ++i__) {
 
260
                                b[i__ + j * b_dim1] += temp * a[i__ + k * 
 
261
                                        a_dim1];
 
262
/* L60: */
 
263
                            }
 
264
                        }
 
265
/* L70: */
 
266
                    }
 
267
/* L80: */
 
268
                }
 
269
            }
 
270
        } else {
 
271
 
 
272
/*           Form  B := alpha*A**T*B. */
 
273
 
 
274
            if (upper) {
 
275
                i__1 = *n;
 
276
                for (j = 1; j <= i__1; ++j) {
 
277
                    for (i__ = *m; i__ >= 1; --i__) {
 
278
                        temp = b[i__ + j * b_dim1];
 
279
                        if (nounit) {
 
280
                            temp *= a[i__ + i__ * a_dim1];
 
281
                        }
 
282
                        i__2 = i__ - 1;
 
283
                        for (k = 1; k <= i__2; ++k) {
 
284
                            temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
 
285
/* L90: */
 
286
                        }
 
287
                        b[i__ + j * b_dim1] = *alpha * temp;
 
288
/* L100: */
 
289
                    }
 
290
/* L110: */
 
291
                }
 
292
            } else {
 
293
                i__1 = *n;
 
294
                for (j = 1; j <= i__1; ++j) {
 
295
                    i__2 = *m;
 
296
                    for (i__ = 1; i__ <= i__2; ++i__) {
 
297
                        temp = b[i__ + j * b_dim1];
 
298
                        if (nounit) {
 
299
                            temp *= a[i__ + i__ * a_dim1];
 
300
                        }
 
301
                        i__3 = *m;
 
302
                        for (k = i__ + 1; k <= i__3; ++k) {
 
303
                            temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
 
304
/* L120: */
 
305
                        }
 
306
                        b[i__ + j * b_dim1] = *alpha * temp;
 
307
/* L130: */
 
308
                    }
 
309
/* L140: */
 
310
                }
 
311
            }
 
312
        }
 
313
    } else {
 
314
        if (igraphlsame_(transa, "N")) {
 
315
 
 
316
/*           Form  B := alpha*B*A. */
 
317
 
 
318
            if (upper) {
 
319
                for (j = *n; j >= 1; --j) {
 
320
                    temp = *alpha;
 
321
                    if (nounit) {
 
322
                        temp *= a[j + j * a_dim1];
 
323
                    }
 
324
                    i__1 = *m;
 
325
                    for (i__ = 1; i__ <= i__1; ++i__) {
 
326
                        b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
 
327
/* L150: */
 
328
                    }
 
329
                    i__1 = j - 1;
 
330
                    for (k = 1; k <= i__1; ++k) {
 
331
                        if (a[k + j * a_dim1] != 0.) {
 
332
                            temp = *alpha * a[k + j * a_dim1];
 
333
                            i__2 = *m;
 
334
                            for (i__ = 1; i__ <= i__2; ++i__) {
 
335
                                b[i__ + j * b_dim1] += temp * b[i__ + k * 
 
336
                                        b_dim1];
 
337
/* L160: */
 
338
                            }
 
339
                        }
 
340
/* L170: */
 
341
                    }
 
342
/* L180: */
 
343
                }
 
344
            } else {
 
345
                i__1 = *n;
 
346
                for (j = 1; j <= i__1; ++j) {
 
347
                    temp = *alpha;
 
348
                    if (nounit) {
 
349
                        temp *= a[j + j * a_dim1];
 
350
                    }
 
351
                    i__2 = *m;
 
352
                    for (i__ = 1; i__ <= i__2; ++i__) {
 
353
                        b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
 
354
/* L190: */
 
355
                    }
 
356
                    i__2 = *n;
 
357
                    for (k = j + 1; k <= i__2; ++k) {
 
358
                        if (a[k + j * a_dim1] != 0.) {
 
359
                            temp = *alpha * a[k + j * a_dim1];
 
360
                            i__3 = *m;
 
361
                            for (i__ = 1; i__ <= i__3; ++i__) {
 
362
                                b[i__ + j * b_dim1] += temp * b[i__ + k * 
 
363
                                        b_dim1];
 
364
/* L200: */
 
365
                            }
 
366
                        }
 
367
/* L210: */
 
368
                    }
 
369
/* L220: */
 
370
                }
 
371
            }
 
372
        } else {
 
373
 
 
374
/*           Form  B := alpha*B*A**T. */
 
375
 
 
376
            if (upper) {
 
377
                i__1 = *n;
 
378
                for (k = 1; k <= i__1; ++k) {
 
379
                    i__2 = k - 1;
 
380
                    for (j = 1; j <= i__2; ++j) {
 
381
                        if (a[j + k * a_dim1] != 0.) {
 
382
                            temp = *alpha * a[j + k * a_dim1];
 
383
                            i__3 = *m;
 
384
                            for (i__ = 1; i__ <= i__3; ++i__) {
 
385
                                b[i__ + j * b_dim1] += temp * b[i__ + k * 
 
386
                                        b_dim1];
 
387
/* L230: */
 
388
                            }
 
389
                        }
 
390
/* L240: */
 
391
                    }
 
392
                    temp = *alpha;
 
393
                    if (nounit) {
 
394
                        temp *= a[k + k * a_dim1];
 
395
                    }
 
396
                    if (temp != 1.) {
 
397
                        i__2 = *m;
 
398
                        for (i__ = 1; i__ <= i__2; ++i__) {
 
399
                            b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
 
400
/* L250: */
 
401
                        }
 
402
                    }
 
403
/* L260: */
 
404
                }
 
405
            } else {
 
406
                for (k = *n; k >= 1; --k) {
 
407
                    i__1 = *n;
 
408
                    for (j = k + 1; j <= i__1; ++j) {
 
409
                        if (a[j + k * a_dim1] != 0.) {
 
410
                            temp = *alpha * a[j + k * a_dim1];
 
411
                            i__2 = *m;
 
412
                            for (i__ = 1; i__ <= i__2; ++i__) {
 
413
                                b[i__ + j * b_dim1] += temp * b[i__ + k * 
 
414
                                        b_dim1];
 
415
/* L270: */
 
416
                            }
 
417
                        }
 
418
/* L280: */
 
419
                    }
 
420
                    temp = *alpha;
 
421
                    if (nounit) {
 
422
                        temp *= a[k + k * a_dim1];
 
423
                    }
 
424
                    if (temp != 1.) {
 
425
                        i__1 = *m;
 
426
                        for (i__ = 1; i__ <= i__1; ++i__) {
 
427
                            b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
 
428
/* L290: */
 
429
                        }
 
430
                    }
 
431
/* L300: */
 
432
                }
 
433
            }
 
434
        }
 
435
    }
 
436
 
 
437
    return 0;
 
438
 
 
439
/*     End of DTRMM . */
 
440
 
 
441
} /* igraphdtrmm_ */
 
442