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 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398
|
c----------------------------------------------------------------------c
subroutine amub (nrow,ncol,job,a,ja,ia,b,jb,ib,
* c,jc,ic,nzmax,iw,ierr)
double precision a(*), b(*), c(*)
integer ja(*),jb(*),jc(*),ia(nrow+1),ib(*),ic(*),iw(ncol)
c-----------------------------------------------------------------------
c performs the matrix by matrix product C = A B
c-----------------------------------------------------------------------
c on entry:
c ---------
c nrow = integer. The row dimension of A = row dimension of C
c ncol = integer. The column dimension of B = column dimension of C
c job = integer. Job indicator. When job = 0, only the structure
c (i.e. the arrays jc, ic) is computed and the
c real values are ignored.
c
c a,
c ja,
c ia = Matrix A in compressed sparse row format.
c
c b,
c jb,
c ib = Matrix B in compressed sparse row format.
c
c nzmax = integer. The length of the arrays c and jc.
c amub will stop if the result matrix C has a number
c of elements that exceeds exceeds nzmax. See ierr.
c
c on return:
c----------
c c,
c jc,
c ic = resulting matrix C in compressed sparse row sparse format.
c
c ierr = integer. serving as error message.
c ierr = 0 means normal return,
c ierr .gt. 0 means that amub stopped while computing the
c i-th row of C with i=ierr, because the number
c of elements in C exceeds nzmax.
c
c work arrays:
c------------
c iw = integer work array of length equal to the number of
c columns in B.
c Note:
c-------
c The row dimension of B is not needed. However there is no checking
c on the condition that ncol(A) = nrow(B).
c
c-----------------------------------------------------------------------
double precision scal
logical values
values = (job .ne. 0)
len = 0
ic(1) = 1
ierr = 0
c initialize array iw.
do 1 j=1, ncol
iw(j) = 0
1 continue
c
do 500 ii=1, nrow
c row i
do 200 ka=ia(ii), ia(ii+1)-1
if (values) scal = a(ka)
jj = ja(ka)
do 100 kb=ib(jj),ib(jj+1)-1
jcol = jb(kb)
jpos = iw(jcol)
if (jpos .eq. 0) then
len = len+1
if (len .gt. nzmax) then
ierr = ii
return
endif
jc(len) = jcol
iw(jcol)= len
if (values) c(len) = scal*b(kb)
else
if (values) c(jpos) = c(jpos) + scal*b(kb)
endif
100 continue
200 continue
do 201 k=ic(ii), len
iw(jc(k)) = 0
201 continue
ic(ii+1) = len+1
500 continue
return
c-------------end-of-amub-----------------------------------------------
c-----------------------------------------------------------------------
end
c-----------------------------------------------------------------------
subroutine amudia (nrow,job, a, ja, ia, diag, b, jb, ib)
double precision a(*), b(*), diag(*)
integer ja(*),jb(*), ia(nrow+1),ib(nrow+1)
c-----------------------------------------------------------------------
c performs the matrix by matrix product B = A * Diag (in place)
c-----------------------------------------------------------------------
c on entry:
c ---------
c nrow = integer. The row dimension of A
c
c job = integer. job indicator. Job=0 means get array b only
c job = 1 means get b, and the integer arrays ib, jb.
c
c a,
c ja,
c ia = Matrix A in compressed sparse row format.
c
c diag = diagonal matrix stored as a vector dig(1:n)
c
c on return:
c----------
c
c b,
c jb,
c ib = resulting matrix B in compressed sparse row sparse format.
c
c Notes:
c-------
c 1) The column dimension of A is not needed.
c 2) algorithm in place (B can take the place of A).
c-----------------------------------------------------------------
do 1 ii=1,nrow
c
c scale each element
c
k1 = ia(ii)
k2 = ia(ii+1)-1
do 2 k=k1, k2
b(k) = a(k)*diag(ja(k))
2 continue
1 continue
c
if (job .eq. 0) return
c
do 3 ii=1, nrow+1
ib(ii) = ia(ii)
3 continue
do 31 k=ia(1), ia(nrow+1) -1
jb(k) = ja(k)
31 continue
return
c-----------------------------------------------------------------------
c-----------end-of-amudiag----------------------------------------------
end
c----------------------------------------------------------------------c
subroutine amux (n, x, y, a,ja,ia)
double precision x(*), y(*), a(*)
integer n, ja(*), ia(*)
c-----------------------------------------------------------------------
c A times a vector
c-----------------------------------------------------------------------
c multiplies a matrix by a vector using the dot product form
c Matrix A is stored in compressed sparse row storage.
c
c on entry:
c----------
c n = row dimension of A
c x = real array of length equal to the column dimension of
c the A matrix.
c a, ja,
c ia = input matrix in compressed sparse row format.
c
c on return:
c-----------
c y = real array of length n, containing the product y=Ax
c
c-----------------------------------------------------------------------
c local variables
c
double precision t
integer i, k
c-----------------------------------------------------------------------
do 100 i = 1,n
c
c compute the inner product of row i with vector x
c
t = 0.0d0
do 99 k=ia(i), ia(i+1)-1
t = t + a(k)*x(ja(k))
99 continue
c
c store result in y(i)
c
y(i) = t
100 continue
c
return
end
c-----------------------------------------------------------------------
subroutine aplb (nrow,ncol,job,a,ja,ia,b,jb,ib,
* c,jc,ic,nzmax,iw,ierr)
double precision a(*), b(*), c(*)
integer ja(*),jb(*),jc(*),ia(nrow+1),ib(nrow+1),ic(nrow+1),
* iw(ncol)
c-----------------------------------------------------------------------
c performs the matrix sum C = A+B.
c-----------------------------------------------------------------------
c on entry:
c ---------
c nrow = integer. The row dimension of A and B
c ncol = integer. The column dimension of A and B.
c job = integer. Job indicator. When job = 0, only the structure
c (i.e. the arrays jc, ic) is computed and the
c real values are ignored.
c
c a,
c ja,
c ia = Matrix A in compressed sparse row format.
c
c b,
c jb,
c ib = Matrix B in compressed sparse row format.
c
c nzmax = integer. The length of the arrays c and jc.
c amub will stop if the result matrix C has a number
c of elements that exceeds exceeds nzmax. See ierr.
c
c on return:
c----------
c c,
c jc,
c ic = resulting matrix C in compressed sparse row sparse format.
c
c ierr = integer. serving as error message.
c ierr = 0 means normal return,
c ierr .gt. 0 means that amub stopped while computing the
c i-th row of C with i=ierr, because the number
c of elements in C exceeds nzmax.
c
c work arrays:
c------------
c iw = integer work array of length equal to the number of
c columns in A.
c
c-----------------------------------------------------------------------
logical values
values = (job .ne. 0)
ierr = 0
len = 0
ic(1) = 1
do 1 j=1, ncol
iw(j) = 0
1 continue
c
do 500 ii=1, nrow
c row i
do 200 ka=ia(ii), ia(ii+1)-1
len = len+1
jcol = ja(ka)
if (len .gt. nzmax) goto 999
jc(len) = jcol
if (values) c(len) = a(ka)
iw(jcol)= len
200 continue
c
do 300 kb=ib(ii),ib(ii+1)-1
jcol = jb(kb)
jpos = iw(jcol)
if (jpos .eq. 0) then
len = len+1
if (len .gt. nzmax) goto 999
jc(len) = jcol
if (values) c(len) = b(kb)
iw(jcol)= len
else
if (values) c(jpos) = c(jpos) + b(kb)
endif
300 continue
do 301 k=ic(ii), len
iw(jc(k)) = 0
301 continue
ic(ii+1) = len+1
500 continue
return
999 ierr = ii
return
c------------end of aplb -----------------------------------------------
c-----------------------------------------------------------------------
end
subroutine csrmsr (n,a,ja,ia,ao,jao,wk,iwk,nnzao,ierr)
double precision a(*),ao(*),wk(n)
integer ia(n+1),ja(*),jao(*),iwk(n+1),nnzao,ierr
c-----------------------------------------------------------------------
c Compressed Sparse Row to Modified - Sparse Row
c Sparse row with separate main diagonal
c-----------------------------------------------------------------------
c converts a general sparse matrix a, ja, ia into
c a compressed matrix using a separated diagonal (referred to as
c the bell-labs format as it is used by bell labs semi conductor
c group. We refer to it here as the modified sparse row format.
c Note: this has been coded in such a way that one can overwrite
c the output matrix onto the input matrix if desired by a call of
c the form
c
c call csrmsr (n, a, ja, ia, a, ja, wk,iwk)
c
c In case ao, jao, are different from a, ja, then one can
c use ao, jao as the work arrays in the calling sequence:
c
c call csrmsr (n, a, ja, ia, ao, jao, ao,jao)
c
c-----------------------------------------------------------------------
c
c on entry :
c---------
c a, ja, ia = matrix in csr format. note that the
c algorithm is in place: ao, jao can be the same
c as a, ja, in which case it will be overwritten on it
c upon return.
c nnzao = the number of non-zero entries in ao, jao
c
c on return :
c-----------
c
c ao, jao = sparse matrix in modified sparse row storage format:
c + ao(1:n) contains the diagonal of the matrix.
c + ao(n+2:nnz) contains the nondiagonal elements of the
c matrix, stored rowwise.
c + jao(n+2:nnz) : their column indices
c + jao(1:n+1) contains the pointer array for the nondiagonal
c elements in ao(n+2:nnz) and jao(n+2:nnz).
c i.e., for i .le. n+1 jao(i) points to beginning of row i
c in arrays ao, jao.
c here nnz = number of nonzero elements+1
c ierr:
c = -1 : length of ao, jao < itpr
c
c work arrays:
c------------
c wk = real work array of length n
c iwk = integer work array of length n+1
c
c notes:
c-------
c Algorithm is in place. i.e. both:
c
c call csrmsr (n, a, ja, ia, ao, jao, ao,jao)
c (in which ao, jao, are different from a, ja)
c and
c call csrmsr (n, a, ja, ia, a, ja, wk,iwk)
c (in which wk, jwk, are different from a, ja)
c are OK.
c--------
c coded by Y. Saad Sep. 1989. Rechecked Feb 27, 1990.
c Modified by Pin Ng on June 11, 2002 to provide warning when
c iptr > nnzao+1
c-----------------------------------------------------------------------
icount = 0
c
c store away diagonal elements and count nonzero diagonal elements.
c
do 1 i=1,n
wk(i) = 0.0d0
iwk(i+1) = ia(i+1)-ia(i)
do 2 k=ia(i),ia(i+1)-1
if (ja(k) .eq. i) then
wk(i) = a(k)
icount = icount + 1
iwk(i+1) = iwk(i+1)-1
endif
2 continue
1 continue
c
c compute total length
c
iptr = n + ia(n+1) - icount
if (iptr .gt. nnzao+1) then
ierr = -1
return
endif
c
c copy backwards (to avoid collisions)
c
do 500 ii=n,1,-1
do 100 k=ia(ii+1)-1,ia(ii),-1
j = ja(k)
if (j .ne. ii) then
ao(iptr) = a(k)
jao(iptr) = j
iptr = iptr-1
endif
100 continue
500 continue
c
c compute pointer values and copy wk(*)
c
jao(1) = n+2
do 600 i=1,n
ao(i) = wk(i)
jao(i+1) = jao(i)+iwk(i+1)
600 continue
return
c------------ end of subroutine csrmsr ---------------------------------
c-----------------------------------------------------------------------
end
|