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
|
subroutine spextr1(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, nel_max, ptr, ierr)
*
* PURPOSE
* perform the sparse matrix extraction B = A(i)
* where i is a vector of indices of length ni
*
* With one vector of index, the extraction operation must be
* meant as if A would be stacked as a big sparse
* column vector...
*
* 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, nel_max, ptr(*), ierr
double precision A_R(*), A_I(*), B_R(*), B_I(*)
* internal vars
integer k, ka, kb, ia, ja, ii
* external functions and subroutines called
integer dicho_search
external dicho_search, sz2ptr
ierr = 0
if ( ni.lt.0 ) then ! the case B = A(:) is not treated by this routine
ierr = -2 ! but via the reshape routine (spmat). So this error
return ! must not appeared
endif
kb = 1
if (A_m .eq. 1) then ! B must be in this case a row vector
do k = 1, ni
ka = dicho_search(i(k), A_icol(1), A_mnel(1) )
if (ka .ne. 0) then
if (kb .gt. nel_max) then ! test memoire
ierr = -1
return
endif
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
B_mnel(1) = kb - 1
else
call sz2ptr(A_mnel, A_m-1, ptr) ! need ptr to acces fastly at the beginning of a row
if (A_n .gt. 1) then
do k = 1, ni
ii = i(k)
ia = mod(ii-1, A_m) + 1
ja = (ii-1)/A_m + 1
ka = dicho_search(ja, A_icol(ptr(ia)), A_mnel(ia) )
if (ka .ne. 0) then
if (kb .gt. nel_max) then ! test memoire
ierr = -1
return
endif
ka = ka + ptr(ia) - 1
B_mnel(k) = 1
B_icol(kb) = 1
if (it .ge. 0) B_R(kb) = A_R(ka)
if (it .eq. 1) B_I(kb) = A_I(ka)
kb = kb + 1
else
B_mnel(k) = 0
endif
enddo
else ! A is a sparse column => binary (dicho) search not useful
do k = 1, ni
ii = i(k)
if ( A_mnel(ii) .gt. 0 ) then
if (kb .gt. nel_max) then ! test memoire
ierr = -1
return
endif
ka = ptr(ii)
B_mnel(k) = 1
B_icol(kb) = 1
if (it .ge. 0) B_R(kb) = A_R(ka)
if (it .eq. 1) B_I(kb) = A_I(ka)
kb = kb + 1
else
B_mnel(k) = 0
endif
enddo
endif
endif
B_nel = kb - 1
end
|
