CoCalc -- mul

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
gap4r8 / pkg / carat / carat-2.1b1 / functions / Matrix / mul_mat.c
⁵⁶³⁶⁸⁰ views
1
#include "typedef.h"
2
#include "tools.h"
3
#include "matrix.h"
4

5
/**************************************************************************\
6
@---------------------------------------------------------------------------
7
@---------------------------------------------------------------------------
8
@ FILE: mul_mat.c
9
@---------------------------------------------------------------------------
10
@---------------------------------------------------------------------------
11
@
12
\**************************************************************************/
13
/*{{{}}}*/
14
/*{{{  rmat_mul, done*/
15
static matrix_TYP *rmat_mul (L_mat, R_mat)
16
matrix_TYP *L_mat, *R_mat ;
17

18
{
19
matrix_TYP *P_mat;
20
int **L, **LN, **R, **RN, **P, **PN ;
21
int  rL, cL, rR, cR, rP, cP ;
22
int i, j, k ; 
23
int tmp_ggt, tmp_z, tmp_n, tmp_kgv;
24

25
  cL =        L_mat->cols;
26
  rL = rP =   L_mat->rows;
27
  rR =        R_mat->rows;
28
  cR = cP =   R_mat->cols;
29
  
30
  P_mat = init_mat(rP,cP,"r");
31
  
32
  L = L_mat->array.SZ;
33
  R = R_mat->array.SZ;
34
  P = P_mat->array.SZ;
35
  LN = L_mat->array.N;
36
  RN = R_mat->array.N;
37
  PN = P_mat->array.N;
38
  
39
  for (i = 0; i < rL; i++)
40
    for (j = 0; j < cL; j++ ) 
41
    {
42
      if ( L[i][j]) {
43
        if ( R_mat->flags.Diagonal ) {
44
          if (R[j][j]) { 
45
            P[i][j]  =  L[i][j] *  R[j][j];
46
            PN[i][j] = LN[i][j] * RN[j][j];
47
          }
48
        } else {
49
          for (k = 0; k < cR; k++) {
50
            if (R[j][k]) { 
51
              tmp_z =  L[i][j] *  R[j][k];
52
              tmp_n = LN[i][j] * RN[j][k];
53
              /*
54
               *   direkt kuerzen! Ueberlaufgefahr
55
               */
56
              Normal2( &tmp_z, &tmp_n );
57
              tmp_ggt = GGT( PN[i][k], tmp_n );
58
              if ( PN[i][k] > tmp_n ) {
59
                tmp_kgv = PN[i][k] / tmp_ggt * tmp_n;
60
              } else {
61
                tmp_kgv = tmp_n/ tmp_ggt * PN[i][k];
62
              }                                   
63
              P[i][k]  = P[i][k] * (tmp_kgv / PN[i][k]);
64
              P[i][k] += tmp_z * (tmp_kgv / tmp_n );
65
              PN[i][k] = tmp_kgv;
66
              Normal2( &P[i][k], &PN[i][k]);
67
            }
68
          }
69
        }           
70
      }
71
    }
72
  
73
  P_mat->kgv = 1;
74
  P_mat->flags.Integral  = FALSE;
75
  P_mat->flags.Scalar = L_mat->flags.Scalar && R_mat->flags.Scalar;
76
  P_mat->flags.Diagonal = L_mat->flags.Diagonal && R_mat->flags.Diagonal;
77
  P_mat->flags.Symmetric = P_mat->flags.Diagonal;
78
  Check_mat( P_mat);
79
  return ( P_mat );
80
}
81
/*}}}  */
82
/*{{{  intmat_mul, done*/
83
static matrix_TYP *intmat_mul (L_mat, R_mat)
84
matrix_TYP *L_mat, *R_mat ;
85

86
{
87
matrix_TYP *P_mat;
88
int **L, **R, **P ;
89
int *L_i,
90
    *R_j,
91
    *P_i;
92
int  rL, cL, rR, cR, rP, cP ;
93
int i, j, k ;
94

95
  cL =        L_mat->cols;
96
  rL = rP =   L_mat->rows;
97
  rR =        R_mat->rows;
98
  cR = cP =   R_mat->cols;
99
  
100
  P_mat = init_mat(rP,cP,"");
101
  
102
  L = L_mat->array.SZ;
103
  R = R_mat->array.SZ;
104
  P = P_mat->array.SZ;
105
  
106
  for (i = 0; i < rL; i++){
107
    L_i = L[i];
108
    P_i = P[i];
109
    for (j = 0; j < cL; j++) 
110
    {
111
      R_j = R[j];
112
      if ( L_i[j])
113
        if ( R_mat->flags.Diagonal ) 
114
        {
115
          if (R_j[j]) P_i[j] = L_i[j] * R_j[j];
116
        }
117
        else
118
          for (k = 0; k < cR; k++)
119
            if (R_j[k])
120
               P_i[k] += L_i[j] * R_j[k];
121
    }
122
  } 
123
  P_mat->kgv = L_mat->kgv * R_mat->kgv;
124
  P_mat->flags.Integral  = (P_mat->kgv == 1);
125
  P_mat->flags.Scalar = L_mat->flags.Scalar && R_mat->flags.Scalar;
126
  P_mat->flags.Diagonal = L_mat->flags.Diagonal && R_mat->flags.Diagonal;
127
  P_mat->flags.Symmetric = P_mat->flags.Diagonal;
128
  Check_mat( P_mat);
129
  return ( P_mat );
130
}
131
/*}}}  */
132
/*{{{  pmat_mul, done*/
133

134

135
/**************************************************************************\
136
@---------------------------------------------------------------------------
137
@ matrix_TYP *pmat_mul (L_mat, R_mat)
138
@ matrix_TYP *L_mat, *R_mat ;
139
@
140
@ calculates L_mat * R_mat modulo L_mat->prime
141
@---------------------------------------------------------------------------
142
@
143
\**************************************************************************/
144
matrix_TYP *pmat_mul (L_mat, R_mat)
145
matrix_TYP *L_mat, *R_mat ;
146
{
147
matrix_TYP *P_mat;
148

149
int **L, rL, cL, cR,  **R, **M ;
150

151
int i, j, k;
152
	
153
rL = L_mat->rows;
154
cL = L_mat->cols;
155
cR = R_mat->cols;
156
P_mat = init_mat(rL,cR,"p");
157
L = L_mat->array.SZ;
158
R = R_mat->array.SZ;
159
M = P_mat->array.SZ;
160
P_mat->prime = L_mat->prime;
161
init_prime(P_mat->prime);
162

163
for ( i = 0 ; i < rL; i++)
164
	for (j = 0; j < cL; j++ )
165
		if ( L[i][j])
166
			for (k = 0 ; k < cR; k++)
167
				if ( R[j][k] )
168
					M[i][k] = S(M[i][k],P(L[i][j],R[j][k]));
169
Check_mat(P_mat);
170
return ( P_mat );
171
}
172

173
/*}}}  */
174
/*{{{  mat_mul, done*/
175

176

177
/**************************************************************************\
178
@---------------------------------------------------------------------------
179
@ matrix_TYP *mat_mul (L_mat, R_mat)
180
@ matrix_TYP *L_mat, *R_mat ;
181
@
182
@ calculates L_mat * R_mat.
183
@---------------------------------------------------------------------------
184
@
185
\**************************************************************************/
186
matrix_TYP *mat_mul (L_mat, R_mat)
187
matrix_TYP *L_mat, *R_mat ;
188

189
{
190
matrix_TYP *P_mat, *tmp_mat;
191

192

193
  if ( L_mat->cols != R_mat->rows ) {
194
    fprintf (stderr, "Can't multiply %dx%d with %dx%d\n", L_mat->rows,
195
                                L_mat->cols, R_mat->rows, R_mat->cols );
196
    exit (3);
197
  }
198
  if ( L_mat->prime != R_mat->prime ) {
199
    fprintf(stderr,"Can`t multiply Fp-matrix and Z-matrix\n");
200
    exit(3);
201
  }
202
  if ( L_mat->prime ) {
203
    P_mat = pmat_mul(L_mat, R_mat);
204
  } else {
205
    if ( L_mat->array.N != NULL ) {
206
      if ( R_mat->array.N != NULL ) {
207
        P_mat = rmat_mul( L_mat, R_mat );
208
      } else {
209
        tmp_mat = copy_mat( R_mat );
210
        kgv2rat( tmp_mat );
211
        P_mat = rmat_mul( L_mat, tmp_mat );
212
        free_mat( tmp_mat );
213
      } 
214
    } else if ( R_mat->array.N != NULL ) {
215
      tmp_mat = copy_mat( L_mat );
216
      kgv2rat( tmp_mat );
217
      P_mat = rmat_mul( tmp_mat, R_mat );
218
      free_mat( tmp_mat );
219
    } else {
220
      P_mat = intmat_mul(L_mat, R_mat);
221
    }
222
  }
223
  return(P_mat);
224
}
225

226
/*}}}  */
227
/*{{{  mat_muleq, unchanged*/
228

229
/**************************************************************************\
230
@---------------------------------------------------------------------------
231
@ matrix_TYP *mat_muleq(L_mat,R_mat)
232
@ matrix_TYP *L_mat, *R_mat;
233
@
234
@ calculates L_mat = L_mat * R_mat.
235
@---------------------------------------------------------------------------
236
@
237
\**************************************************************************/
238
matrix_TYP *mat_muleq(L_mat,R_mat)
239
matrix_TYP *L_mat, *R_mat;
240

241
{
242
matrix_TYP *H_mat;
243
matrix_TYP merk;
244

245
  H_mat = mat_mul(L_mat,R_mat);
246

247
  memcpy( &merk, L_mat, sizeof(matrix_TYP) );
248
  memcpy( L_mat, H_mat, sizeof(matrix_TYP) );
249

250
  /* inserted 21/1/97 tilman */
251
  memcpy( H_mat, &merk, sizeof(matrix_TYP) );
252

253

254
  free_mat(H_mat);
255

256
  return(L_mat);
257
}
258
/*}}}  */
259
/*{{{  mat_kon*/
260

261

262
/**************************************************************************\
263
@---------------------------------------------------------------------------
264
@ matrix_TYP *mat_kon(L_mat,M_mat,R_mat)
265
@ matrix_TYP *L_mat, *M_mat, *R_mat;
266
@
267
@ calculates L_mat * M_mat * R_mat
268
@---------------------------------------------------------------------------
269
@
270
\**************************************************************************/
271
matrix_TYP *mat_kon(L_mat,M_mat,R_mat)
272
matrix_TYP *L_mat, *M_mat, *R_mat;
273
{
274
matrix_TYP *P_mat, *T_mat;
275

276

277
  if ( L_mat->cols != M_mat->rows ) 
278
  {
279
    fprintf (stderr, "Can't multiply %dx%d with %dx%d\n", L_mat->rows,
280
                                L_mat->cols, M_mat->rows, M_mat->cols );
281
    exit (3);
282
  }
283
  if ( M_mat->cols != R_mat->rows ) 
284
  {
285
    fprintf (stderr, "Can't multiply %dx%d with %dx%d\n", L_mat->rows,
286
                                L_mat->cols, R_mat->rows, R_mat->cols );
287
    exit (3);
288
  }
289
  if ( L_mat->prime != R_mat->prime )
290
  {
291
    fprintf(stderr,"Can`t multiply Fp-matrix and Z-matrix\n");
292
    exit(3);
293
  }
294
  if ( L_mat->prime != M_mat->prime )
295
  {
296
    fprintf(stderr,"Can`t multiply Fp-matrix and Z-matrix\n");
297
    exit(3);
298
  }
299
  T_mat = mat_mul(L_mat, M_mat);
300
  P_mat = mat_mul(T_mat, R_mat);
301
  free_mat(T_mat);
302
  return(P_mat);
303
}
304

305
/*}}}  */
306

307
Product

Resources

Company
