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
|
subroutine spif(A_m, A_n, A_nel, A_it, A_mnel, A_icol, A_R, A_I,
$ B_m, B_n, B_it, B_R, B_I,
$ C_m, C_n, C_nel, C_it, C_mnel, C_icol, C_R, C_I,
$ i, pi, ni, j, pj, nj, nelmax, ierr)
*
* PURPOSE
* do the scilab operation A(i,j) = B : where A is a sparse matrix
* (real if it = 0, complex if it=1), B is a full matrix and i and
* j are 2 arrays of indices i(1..ni), j(1..nj).
* The result being stored in C (a sparse). nelmax is the max number of
* elements for the resulting matrix C.
*
* AUTHOR
* Bruno Pincon
*
implicit none
integer A_m, A_n, A_nel, A_it, A_mnel(*), A_icol(*),
$ B_m, B_n, B_it,
$ C_m, C_n, C_nel, C_it, C_mnel(*), C_icol(*),
$ i(*), pi(*), ni, j(*), pj(*), nj, nelmax, ierr
double precision A_R(*), A_I(*), B_R(B_m,*), B_I(B_m,*),
$ C_R(*), C_I(*)
* local var
logical allrow, allcol, B_is_scalar
integer ii, i1, i2, imax, k, ka, kc, A_mnel_ii
* external functions and subroutines
logical is_in_order
external is_in_order, isorti, set_perm_id, copy_fullrow2sprow,
$ insert_row, iset, copy_sprow
ierr = 0
allrow = ni.lt.0 ! ni < 0 for A(:,j) = B
allcol = nj.lt.0 ! nj < 0 for A(i,:) = B
if (allrow) then
ni = A_m
else
if ( .not. is_in_order(i, ni) ) then
call isorti(i, pi, ni)
else
call set_perm_id(pi, ni)
endif
imax = i(pi(ni))
endif
if (allcol) then
nj = A_n
else
if ( .not. is_in_order(j, nj) ) then
call isorti(j, pj, nj)
else
call set_perm_id(pj, nj)
endif
endif
if ( (B_m .eq. 1) .and. (B_n .eq.1) ) then
B_is_scalar = .true.
B_m = ni
B_n = nj
else
B_is_scalar = .false.
endif
ka = 1
kc = 1
if (allrow .and. allcol) then ! *** A(:,:) = B
do ii = 1, B_m
call copy_fullrow2sprow(ii, kc, C_it, C_mnel(ii), C_icol,
$ C_R, C_I, B_m, B_n, B_it, B_R, B_I,
$ B_is_scalar, nelmax, ierr)
if (ierr .ne. 0) return
enddo
elseif (allrow) then ! *** A(:,j) = B
call iset(A_m, 0, C_mnel, 1)
do ii = 1, A_m
call insert_row(ka, A_it, A_mnel(ii), A_icol, A_R, A_I,
$ kc, C_it, C_mnel(ii), C_icol, C_R, C_I,
$ j, pj, nj, ii, B_m, B_it, B_R, B_I,
$ B_is_scalar, nelmax, ierr)
if (ierr .ne. 0) return
enddo
else ! *** A(i,:) = B or A(i,j) = B
call iset(imax, 0, C_mnel, 1)
i1 = 1
k = 1
* loop on k from 1 to ni but some values must be skiped
100 ii = i(pi(k))
if (k .lt. ni) then
if (i(pi(k+1)) .eq. ii) then
k = k+1
goto 100
endif
endif
i2 = min(ii-1, A_m)
call copy_sprow(i1,i2, ka, A_it, A_mnel, A_icol, A_R, A_I,
$ kc, C_it, C_mnel, C_icol, C_R, C_I,
$ nelmax, ierr)
if (ierr .ne. 0) return
if (allcol) then ! *** A(i,:) = B
call copy_fullrow2sprow(pi(k), kc, C_it, C_mnel(ii),
$ C_icol, C_R, C_I,
$ B_m, B_n, B_it, B_R, B_I,
$ B_is_scalar, nelmax, ierr)
ka = ka + A_mnel(ii)
else ! *** A(i,j) = B
if ( ii .gt. A_m ) then
A_mnel_ii = 0
else
A_mnel_ii = A_mnel(ii)
endif
call insert_row(ka, A_it, A_mnel_ii, A_icol, A_R, A_I,
$ kc, C_it, C_mnel(ii), C_icol, C_R, C_I,
$ j, pj, nj, pi(k), B_m, B_it, B_R, B_I,
$ B_is_scalar, nelmax, ierr)
endif
if (ierr .ne. 0) return
i1 = ii+1
k = k+1
if (k .le. ni) goto 100
* endloop on k
i1 = i(pi(ni))+1
i2 = A_m
call copy_sprow(i1,i2, ka, A_it, A_mnel, A_icol, A_R, A_I,
$ kc, C_it, C_mnel, C_icol, C_R, C_I,
$ nelmax, ierr)
endif
C_nel = kc-1
end
|
