1from paida.paida_core.PAbsorber import * 2from paida.math.pylapack.pyblas.xerbla import xerbla 3from paida.math.pylapack.pyblas.lsame import lsame 4 5def dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc): 6 """Written on 8-February-1989. 7 Jack Dongarra, Argonne National Laboratory. 8 Iain Duff, AERE Harwell. 9 Jeremy Du Croz, Numerical Algorithms Group Ltd. 10 Sven Hammarling, Numerical Algorithms Group Ltd. 11 Python replacement by K. KISHIMOTO (korry@users.sourceforge.net) 12 13 Level 3 Blas routine. 14 15 Purpose 16 ======= 17 18 dgemm performs one of the matrix-matrix operations 19 C := alpha*op( A )*op( B ) + beta*C, 20 where op( X ) is one of 21 op( X ) = X or op( X ) = X', 22 alpha and beta are scalars, 23 and A, B and C are matrices, 24 with op( A ) an m by k matrix, 25 op( B ) a k by n matrix 26 and C an m by n matrix. 27 28 Parameters 29 ========== 30 31 transa - character 32 On entry, transa specifies the form of op( A ) 33 to be used in the matrix multiplication as follows: 34 transa = 'N' or 'n', op( A ) = A. 35 transa = 'T' or 't', op( A ) = A'. 36 transa = 'C' or 'c', op( A ) = A'. 37 Unchanged on exit. 38 39 transb - character 40 On entry, transb specifies the form of op( B ) 41 to be used in the matrix multiplication as follows: 42 transb = 'N' or 'n', op( B ) = B. 43 transb = 'T' or 't', op( B ) = B'. 44 transb = 'C' or 'c', op( B ) = B'. 45 Unchanged on exit. 46 47 m - integer 48 On entry, m specifies the number of rows of the matrix op( A ) and of the matrix C. 49 m must be at least zero. 50 Unchanged on exit. 51 52 n - integer 53 On entry, n specifies the number of columns of the matrix op( B ) 54 and of the number of columns of the matrix C. 55 n must be at least zero. 56 Unchanged on exit. 57 58 k - integer 59 On entry, k specifies the number of columns of the matrix op( A ) 60 and the number of rows of the matrix op( B ). 61 k must be at least zero. 62 Unchanged on exit. 63 64 alpha - double 65 On entry, alpha specifies the scalar alpha. 66 Unchanged on exit. 67 68 a - double array of dimension at least (min(k, m), min(k, m)). 69 Before entry with transa[0] = 'N' or 'n', 70 the leading m by k part of the array A must contain the matrix A, 71 otherwise the leading k by m part of the array A must contain the matrix A. 72 Unchanged on exit. 73 74 lda - leading dimension. 75 In almost all cases, it is for compatibility. 76 77 b - double array of dimension at least (min(k, n), min(k, n)). 78 Before entry with transa[0] = 'N' or 'n', 79 the leading n by k part of the array b must contain the matrix B, 80 otherwise the leading k by n part of the array b must contain the matrix B. 81 Unchanged on exit. 82 83 ldb - leading dimension. 84 In almost all cases, it is for compatibility. 85 86 beta - double 87 On entry, beta specifies the scalar beta. 88 When beta is supplied as zero then C need not be set on input. 89 Unchanged on exit. 90 91 C - double array of dimension at least (max(l, m), n). 92 Before entry, the leading m by n part of the array C must contain the matrix C, 93 except when beta is zero, in which case C need not be set on entry. 94 On exit, the array C is overwritten by the m by n matrix ( alpha*op( A )*op( B ) + beta*C ). 95 96 ldc - leading dimension. 97 In almost all cases, it is for compatibility. 98 99 External Subroutines 100 ==================== 101 102 lsame, xerbla 103 """ 104 105 one = 1.0e+0 106 zero = 0.0e+0 107 108 ### Set NOTA and NOTB as true if A and B respectively are not transposed 109 ### and set NROWA, NCOLA and NROWB as the number of rows 110 ### and columns of A and the number of rows of B respectively. 111 nota = lsame(transa, 'N') 112 notb = lsame(transb, 'N') 113 if nota: 114 nrowa = m 115 ncola = k 116 else: 117 nrowa = k 118 ncola = m 119 if notb: 120 nrowb = k 121 else: 122 nrowb = n 123 124 ### Test the input parameters. 125 info[0] = 0 126 if (not nota) and (not lsame(transa, 'C')) and (not lsame(transa, 'T')): 127 info[0] = 1 128 elif (not notb) and (not lsame(transb, 'C')) and (not lsame(transb, 'T')): 129 info[0] = 2 130 elif m < 0: 131 info[0] = 3 132 elif n < 0: 133 info[0] = 4 134 elif k < 0: 135 info[0] = 5 136 if info[0] != 0: 137 xerbla('dgemm', info) 138 return 139 140 ### Quick return if possible. 141 if (m == 0) or (n == 0) or (((alpha == 0) or (k == 0)) and (beta == one)): 142 return 143 144 ### And if alpha == zero. 145 if alpha == zero: 146 if beta == zero: 147 for j in range(1, n + 1): 148 for i in range(1, m + 1): 149 c[j - 1, i - 1] = zero 150 else: 151 for j in range(1, n + 1): 152 for i in range(1, m + 1): 153 c[j - 1, i - 1] *= beta 154 return 155 156 ### Start the operations. 157 if notb: 158 if nota: 159 ### Form C := alpha*A*B + beta*C. 160 for j in range(1, n + 1): 161 if beta == zero: 162 for i in range(1, m + 1): 163 c[j - 1, i - 1] = zero 164 elif beta != one: 165 for i in range(1, m + 1): 166 c[j - 1, i - 1] *= beta 167 for l in range(1, k + 1): 168 if b[j - 1, l - 1] != zero: 169 temp = alpha * b[j - 1, l - 1] 170 for i in range(1, m + 1): 171 c[j - 1, i - 1] += temp * a[l - 1, i - 1] 172 else: 173 ### Form C := alpha*A'*B + beta*C 174 for j in range(1, n + 1): 175 for i in range(1, m + 1): 176 temp = zero 177 for l in range(1, k + 1): 178 temp += a[i - 1, l - 1] * b[j - 1, l - 1] 179 if beta == zero: 180 c[j - 1, i - 1] = alpha * temp 181 else: 182 c[j - 1, i - 1] = alpha * temp + beta * c[j - 1, i - 1] 183 else: 184 if nota: 185 ### Form C := alpha*A*B' + beta*C 186 for j in range(1, n + 1): 187 if beta == zero: 188 for i in range(1, m + 1): 189 c[j - 1, i - 1] = zero 190 elif beta != one: 191 for i in range(1, m + 1): 192 c[j - 1, i - 1] *= beta 193 for l in range(1, k + 1): 194 if b[l - 1, j - 1] != zero: 195 temp = alpha * b[l - 1, j - 1] 196 for i in range(1, m + 1): 197 c[j - 1, i - 1] += temp * a[l - 1, i - 1] 198 else: 199 ### Form C := alpha*A'*B' + beta*C 200 for j in range(1, n + 1): 201 for i in range(1, m + 1): 202 temp = zero 203 for l in range(1, k + 1): 204 temp += a[i - 1, l - 1] * b[l - 1, j - 1] 205 if beta == zero: 206 c[j - 1, i - 1] = alpha * temp 207 else: 208 c[j - 1, i - 1] = alpha * temp + beta * c[j - 1, i - 1] 209 return 210