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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
|
subroutine wspdsp(ne,ind,xr,xi,m,n,maxc,mode,ll,lunit,cw)
C !but
C wspdsp visualise une matrice creuse complexe
C !liste d'appel
C
C subroutine wspdsp(ne,ind,xr,xi,m,n,maxc,mode,ll,lunit,cw)
C
C double precision xr(ne),xi(ne)
C integer ind(*)
C integer nx,m,n,maxc,mode,ll,lunit
C character cw*(*)
C
C c : nombre d'elements nons nuls de la matrice
C ind : indices specifiant la position des elements non nuls
C xr,xi : tableaux contenant les parties reelles et imaginaires des
C elements non nuls
C m : nombre de ligne de la matrice
C n : nombre de colonnes de la matrice
C maxc : nombre de caracteres maximum autorise pour
C representer un nombre
C mode : si mode=1 representation variable
C si mode=0 representation d(maxc).(maxc-7)
C ll : longueur de ligne maximum admissible
C lunit : etiquette logique du support d'edition
C cw : chaine de caracteres de travail de longueur au moins ll
C !
c Copyright INRIA
double precision xr(ne),xi(ne),a,a1,a2,fact,eps,dlamch,ar,ai
integer maxc,mode,fl,typ
integer ind(*)
character cw*(*),sgn*1,dl*1
character*10 form(2)
C
if (ne .eq. 0) then
write (cw,'(''('',i5,'','',i5,'') zero sparse matrix'')') m, n
call basout(io,lunit,cw(1:32))
call basout(io,lunit,' ')
goto 99
else
write (cw,'(''('',i5,'','',i5,'') sparse matrix'')') m, n
call basout(io,lunit,cw(1:27))
call basout(io,lunit,' ')
if (io .eq. -1) goto 99
endif
ilr = 1
ilc = m + 1
nx = 1
eps = dlamch('p')
cw = ' '
write (form(1),130) maxc, maxc-7
dl = ' '
if (m*n .gt. 1) dl = ' '
C
C facteur d'echelle
C
fact = 1.0d+0
a1 = 0.0d+0
if (ne .eq. 1) goto 10
a2 = abs(xr(1)) + abs(xi(1))
do 5 i = 1,ne
a = abs(xr(i)) + abs(xi(i))
if (a.eq.0.0d+0 .or. a.gt.dlamch('o')) goto 5
a1 = max(a1,a)
a2 = min(a2,a)
5 continue
imax = 0
imin = 0
if (a1 .gt. 0.0d+0) imax = int(log10(a1))
if (a2 .gt. 0.0d+0) imin = int(log10(a2))
if (imax*imin .le. 0) goto 10
imax = (imax+imin) / 2
if (abs(imax) .ge. maxc-2) fact = 10.0d+0**(-imax)
10 continue
eps = a1 * fact * eps
C
if (fact .ne. 1.0d+0) then
write (cw(1:12),'(1x,1pd9.1,'' *'')') 1.0d+0/fact
call basout(io,lunit,cw(1:12))
call basout(io,lunit,' ')
if (io .eq. -1) goto 99
endif
i0 = 0
i1 = i0
l = 1
do 20 k = 1,ne
cw = ' '
11 i0 = i0 + 1
if (i0-i1 .gt. ind(l)) then
i1 = i0
l = l + 1
goto 11
endif
i = l
j = ind(ilc-1+k)
write (cw,'(''('',i5,'','',i5,'')'')') i, j
l1 = 18
ar = xr(k) * fact
ai = xi(k) * fact
c if (abs(ar).lt.eps .and. mode.ne.0) ar = 0.0d+0
c if (abs(ai).lt.eps .and. mode.ne.0) ai = 0.0d+0
if (ar.eq.0.0d0 .and. ai.ne.0.0d0) goto 16
C determination du format devant representer a
typ = 1
if (mode .eq. 1) call fmt(abs(ar),maxc,typ,n1,n2)
if (typ .eq. 2) then
fl = n1
ifmt = n2 + 32*n1
elseif (typ .lt. 0) then
ifmt = typ
fl = 3
else
ifmt = 1
fl = maxc
n2 = maxc - 7
endif
sgn = ' '
if (ar .lt. 0.0d+0) sgn = '-'
ar = abs(ar)
cw(l1:l1+1) = ' ' // sgn
l1 = l1 + 2
if (ifmt .eq. 1) then
nf = 1
fl = maxc
n2 = 1
write (cw(l1:l1+fl-1),form(nf)) ar
elseif (ifmt .ge. 0) then
nf = 2
n1 = ifmt / 32
n2 = ifmt - 32*n1
fl = n1
write (form(nf),120) fl, n2
write (cw(l1:l1+fl-1),form(nf)) ar
elseif (ifmt .eq. -1) then
C Inf
fl = 3
cw(l1:l1+fl-1) = 'Inf'
elseif (ifmt .eq. -2) then
C Nan
fl = 3
cw(l1:l1+fl-1) = 'Nan'
endif
l1 = l1 + fl
if (ll .eq. 2) then
cw(l1:l1) = 'i'
l1 = l1 + 1
endif
16 continue
if (ai .eq. 0.0d0) goto 17
typ = 1
if (mode .eq. 1) call fmt(abs(ai),maxc,typ,n1,n2)
if (typ .eq. 2) then
fl = n1
ifmt = n2 + 32*n1
elseif (typ .lt. 0) then
ifmt = typ
fl = 3
else
ifmt = 1
fl = maxc
n2 = maxc - 7
endif
if (ar .ne. 0.0d+0) then
sgn = '+'
else
sgn = ' '
endif
if (ai .lt. 0.0d+0) sgn = '-'
ai = abs(ai)
cw(l1:l1+1) = ' ' // sgn
l1 = l1 + 2
if (ifmt .eq. 1) then
nf = 1
fl = maxc
n2 = 1
write (cw(l1:l1+fl-1),form(nf)) ai
elseif (ifmt .ge. 0) then
nf = 2
n1 = ifmt / 32
n2 = ifmt - 32*n1
fl = n1
write (form(nf),120) fl, n2
write (cw(l1:l1+fl-1),form(nf)) ai
elseif (ifmt .eq. -1) then
C Inf
fl = 3
cw(l1:l1+fl-1) = 'Inf'
elseif (ifmt .eq. -2) then
C Nan
fl = 3
cw(l1:l1+fl-1) = 'Nan'
endif
l1 = l1 + fl
cw(l1:l1) = 'i'
l1 = l1 + 1
17 continue
call basout(io,lunit,cw(1:l1))
if (io .eq. -1) goto 99
20 continue
99 continue
C
120 format ('(f',i2,'.',i2,')')
130 format ('(1pd',i2,'.',i2,')')
end
|
