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
|
subroutine spextr(A_m, A_n, A_nel, A_mnel, A_icol, A_R, A_I,
$ B_m, B_n, B_nel, B_mnel, B_icol, B_R, B_I,
$ it, i, ni, j, nj, nel_max, ptr, p, ierr)
*
* PURPOSE
* perform a sparse matrix extraction, ie, computes the
* sparse matrix B such that : B = A(i,j) where i and
* j are vectors of indices of length ni and nj respectively.
*
* depending upon the argument it, A may be (and so B) :
* a/ a boolean sparse (it = -1). In this case the arrays
* A_R, A_I, B_R, B_I are not used.
* b/ a real sparse (it = 0). The arrays A_I and B_I are not used.
*
* c/ a complex sparse (it = 1)
*
* CAUTION
* nel_max is the maximum non zeros elements authorized
* for B => if this limit is not enought the error
* indicator ierr is set to -1 (else 0)
*
* AUTHOR
* Bruno Pincon
*
implicit none
integer A_m, A_n, A_nel, A_mnel(*), A_icol(*)
integer B_m, B_n, B_nel, B_mnel(*), B_icol(*)
integer it, i(*), ni, j(*), nj, nel_max, ptr(*), p(*), ierr
double precision A_R(*), A_I(*), B_R(*), B_I(*)
* internal vars
integer k, kk, ka, kb, ib, ia, ja, ptrb
logical allrow, allcol, j_in_order
* external functions and subroutines called
integer dicho_search, dicho_search_bis
logical is_in_order
external dicho_search, dicho_search_bis, is_in_order, qsorti,
$ insert_in_order, icopy, unsfdcopy, sz2ptr
ierr = 0
allrow = ni.lt.0 ! ni < 0 for a A(:,j) extraction
allcol = nj.lt.0 ! nj < 0 for a A(i,:) extraction
if (allrow) then
ni = A_m
endif
if (allcol) then
nj = A_n
else
j_in_order = is_in_order(j, nj)
if ( .not. j_in_order ) call qsorti(j, p, nj)
endif
B_m = ni
B_n = nj
call sz2ptr(A_mnel, A_m, ptr)
kb = 1
do ib = 1, B_m
B_mnel(ib) = 0
if (allrow) then
ia = ib
else
ia = i(ib)
endif
if ( A_mnel(ia) .gt. 0 ) then
if (allcol) then
if ( kb+A_mnel(ia) .ge. nel_max ) then ! not enought memory
ierr = -1
return
endif
B_mnel(ib) = A_mnel(ia)
call icopy(A_mnel(ia), A_icol(ptr(ia)), 1, B_icol(kb), 1)
if ( it .ge. 0 )
$ call unsfdcopy(A_mnel(ia), A_R(ptr(ia)), 1, B_R(kb), 1)
if ( it .eq. 1 )
$ call unsfdcopy(A_mnel(ia), A_I(ptr(ia)), 1, B_I(kb), 1)
kb = kb + A_mnel(ia)
elseif ( nj .gt. A_mnel(ia) .and. j_in_order ) then
! algorithm with loop on the row ia of A
do ka = ptr(ia), ptr(ia+1)-1
ja = A_icol(ka)
k = dicho_search(ja, j, nj)
if ( k .ne. 0 ) then ! we have found (the smallest) k such that ja = j(k)
100 if (kb .gt. nel_max) then ! not enought memory
ierr = -1
return
endif
B_mnel(ib) = B_mnel(ib) + 1
B_icol(kb) = k
if ( it .ge. 0 ) B_R(kb) = A_R(ka)
if ( it .eq. 1 ) B_I(kb) = A_I(ka)
kb = kb + 1
if ( k .lt. nj) then ! verify if j(k+1) = ja
k = k + 1
if ( j(k) .eq. ja ) goto 100
endif
endif
enddo
elseif ( nj .gt. 2*A_mnel(ia) .and. .not.j_in_order ) then
! algorithm with loop on the row ia of A but the array
! j is not in order : nevertheless we may work with j(p(.)) where
! p is the permutation which reorders the array j
ptrb = kb
do ka = ptr(ia), ptr(ia+1)-1
ja = A_icol(ka)
k = dicho_search_bis(ja, j, p, nj)
if ( k .ne. 0 ) then ! we have found (the smallest) k such that ja = j(p(k))
200 if (kb .gt. nel_max) then ! not enought memory
ierr = -1
return
endif
B_mnel(ib) = B_mnel(ib) + 1
call insert_in_order( B_icol, ptrb, kb, p(k), it,
$ B_R, B_I, A_R(ka), A_I(ka) )
kb = kb + 1
if ( k .lt. nj) then ! verify if j(p(k+1)) = ja
k = k + 1
if ( j(p(k)) .eq. ja ) goto 200
endif
endif
enddo
else
! algorithm with loop on j(1),...,j(nj)
do k = 1, nj
kk = dicho_search(j(k), A_icol(ptr(ia)), A_mnel(ia))
if ( kk .ne. 0 ) then ! we have found kk such that j(k) = A_icol(ptr(ia)-1+kk)
ka = kk + ptr(ia) - 1
if (kb .gt. nel_max) then ! not enought memory
ierr = -1
return
endif
B_mnel(ib) = B_mnel(ib) + 1
B_icol(kb) = k
if ( it .ge. 0 ) B_R(kb) = A_R(ka)
if ( it .eq. 1 ) B_I(kb) = A_I(ka)
kb = kb + 1
endif
enddo
endif
endif
enddo
B_nel = kb - 1
end
|
