1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
c Copyright (C) 1985-2008 - INRIA - Carlos KLIMANN
c
c This file must be used under the terms of the CeCILL.
c This source file is licensed as described in the file COPYING, which
c you should have received as part of this distribution. The terms
c are also available at
c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
C/MEMBR ADD NAME=DMPMU,SSI=0
c Copyright INRIA
subroutine dmpmu(mp1,d1,nl1,mp2,d2,nl2,mp3,d3,l,m,n)
c!purpose
c ce sous programme effectue le calcul du produit de matrices
c de polynomes a coefficients reels
c
c mp3 = mp1 * mp2
c
c! calling sequence
c
c mp1 : tableau contenant les coefficients des polynomes,
c le coefficient de degre k du polynome mp1(i,j) est range
c dans mp1( d1(i + (j-1)*nl1 + k) )
c mp1 doit etre de taille au moins d1(nl1*n+1)-d1(1)
c d1 : tableau entier de taille nl1*n+1, si k=i+(j-1)*nl1 alors
c d1(k)) contient l'adresse dans mp1 du coeff de degre 0
c du polynome mp1(i,j). Le degre du polynome mp1(i,j) vaut:
c d1(k+1)-d1(k) -1
c nl1 : entier definissant le rangement dans d1
c
c mp2,d2,nl2 : definitions similaires a celles de mp1,d1,nl1
c mp3,d3 : definitions similaires a celles de mp1 et d1, nl3 est
c suppose egal a l
c l: nombre de lignes de la matrice pm1
c m : nombre de ligne des matrices pm
c n : nombre de colonnes des matrices pm
c avec les conventions suivantes:
c
c -si l = 0 signifie que mp1 est un polynome, et que la
c multiplication s'entend alors au sens de la multiplication
c de tous les coefficients de mp2 par mp1.
c
c -si n = 0 signifie que mp2 est un polynome, et que la
c multiplication s'entend au sens de la multiplication de tous
c les coefficients de mp1 par mp2.
c
c -si m = 0 la multiplication s'entend element par element,
c il est alors suppose que mp1 et mp2, donc mp3 sont de dimension
c l x n.
c
c!
c
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
c
c
double precision mp1(*), mp2(*), mp3(*)
integer d1(*), d2(*), d3(*)
integer l, m, n
integer nl1, nl2
integer i, j, k, k1, k2, k3
integer p1, p2, p3
c
d3(1)= 1
if (l .eq. 0 .or. m .eq. 0 .or. n .eq. 0) goto 500
c
p2 = -nl2
p3 = -l
do 10 j = 1, n
p2 = p2 + nl2
p3 = p3 + l
do 11 i = 1, l
mp3(d3(p3+i)) = 0.0d+0
k3 = 0
p1 = i - nl1
do 12 k = 1, m
p1 = p1 + nl1
k2 = d2(p2+k+1) - d2(p2+k) - 1
k1 = d1(p1 + 1) - d1(p1) - 1
call dpmul(mp1(d1(p1)),k1,mp2(d2(p2+k)),k2,
& mp3(d3(p3+i)),k3)
12 continue
d3(p3 + i + 1) = d3(p3 + i) + k3 + 1
11 continue
10 continue
return
500 if (l .eq. 0) goto 600
if (m .eq. 0) goto 700
p1 = -nl1
p3 = -l
k2 = d2(2) - d2(1) - 1
do 510 j = 1, m
p1 = p1 + nl1
p3 = p3 + l
do 510 i = 1, l
k3 = 0
k1 = d1(p1 + i + 1) - d1(p1 + i) - 1
mp3(d3(p3 + i)) = 0.0d+0
call dpmul(mp1(d1(p1+i)),k1,mp2(1),k2,mp3(d3(p3+i)),k3)
d3(p3 + i + 1) = d3(p3 + i) + k3 + 1
510 continue
return
600 k1 = d1(2) - d1(1) - 1
p2 = -nl2
p3 = -m
do 610 j = 1, n
p2 = p2 + nl2
p3 = p3 + m
do 610 i = 1, m
k3 = 0
k2 = d2(p2 + i + 1) - d2(p2 + i) - 1
mp3(d3(p3+i)) = 0.0d+0
call dpmul(mp1(1),k1,mp2(d2(p2+i)),k2,mp3(d3(p3+i)),k3)
d3(p3 + i + 1) = d3(p3 + i) + k3 + 1
610 continue
return
700 p1 = -nl1
p2 = -nl2
p3 = -l
do 710 j = 1, n
p1 = p1 + nl1
p2 = p2 + nl2
p3 = p3 + l
do 710 i = 1, l
k1 = d1(p1 + i + 1) - d1(p1 + i) - 1
k2 = d2(p2 + i + 1) - d2(p2 + i) - 1
mp3(d3(p3+i)) = 0.0d+0
k3 = 0
call dpmul(mp1(d1(p1+i)),k1,mp2(d2(p2+i)),k2,mp3(d3(p3+i)),
& k3)
d3(p3 + i + 1) = d3(p3 + i) + k3 + 1
710 continue
return
end
|
