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
|
program markov
c-----------------------------------------------------------------------
c
c program to generate a Markov chain matrix (to test eigenvalue routines
c or algorithms for singular systems (in which case use I-A ))
c the matrix models simple random walk on a triangular grid.
c see additional comments in subroutine.
c -----
c just compile this segment and link to the rest of sparskit
c (uses subroutine prtmt from MATGEN)
c will create a matrix in the HARWELL/BOEING format and put it in
c the file markov.mat
c
c-----------------------------------------------------------------------
parameter (nmax=5000, nzmax= 4*nmax)
real*8 a(nzmax)
integer ja(nzmax), ia(nmax+1)
c
character title*72,key*8,type*3
open (unit=11,file='markov.mat')
c
c read - in grid size - will not accept too large grids.
c
write (6,'(17hEnter grid-size: ,$)')
read *, m
if (m*(m+1) .gt. 2*nmax ) then
print *, ' m too large - unable to produce matrix '
stop
endif
c
c call generator.
c
call markgen (m, n, a, ja, ia)
c-----------------------------------------------------------------------
title=' Test matrix from SPARSKIT - markov chain model '
key = 'randwk01'
type = 'rua'
iout = 11
job = 2
ifmt = 10
call prtmt (n, n, a, ja, ia, x,'NN',title,
* key, type,ifmt, job, iout)
stop
end
c
subroutine markgen (m, n, a, ja, ia)
c-----------------------------------------------------------------------
c matrix generator for a markov model of a random walk on a triang. grid
c-----------------------------------------------------------------------
c this subroutine generates a test matrix that models a random
c walk on a triangular grid. This test example was used by
c G. W. Stewart ["{SRRIT} - a FORTRAN subroutine to calculate the
c dominant invariant subspaces of a real matrix",
c Tech. report. TR-514, University of Maryland (1978).] and in a few
c papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34,
c pp. 269-295 (1980) ]. These matrices provide reasonably easy
c test problems for eigenvalue algorithms. The transpose of the
c matrix is stochastic and so it is known that one is an exact
c eigenvalue. One seeks the eigenvector of the transpose associated
c with the eigenvalue unity. The problem is to calculate the
c steady state probability distribution of the system, which is
c the eigevector associated with the eigenvalue one and scaled in
c such a way that the sum all the components is equal to one.
c-----------------------------------------------------------------------
c parameters
c------------
c on entry :
c----------
c m = integer. number of points in each direction.
c
c on return:
c----------
c n = integer. The dimension of the matrix. (In fact n is known
c to be equal to (m(m+1))/2 )
c a,
c ja,
c ia = the matrix stored in CSR format.
c
c-----------------------------------------------------------------------
c Notes: 1) the code will actually compute the transpose of the
c stochastic matrix that contains the transition probibilities.
c 2) It should also be possible to have a matrix generator
c with an additional parameter (basically redefining `half' below
c to be another parameter and changing the rest accordingly, but
c this is not as simple as it sounds). This is not likely to provide
c any more interesting matrices.
c-----------------------------------------------------------------------
real*8 a(*), cst, pd, pu, half
integer ja(*), ia(*)
c-----------------------------------------------------------------------
data half/0.5d0/
c
cst = half/real(m-1)
c
c --- ix counts the grid point (natural ordering used), i.e.,
c --- the row number of the matrix.
c
ix = 0
jax = 1
ia(1) = jax
c
c sweep y coordinates
c
do 20 i=1,m
jmax = m-i+1
c
c sweep x coordinates
c
do 10 j=1,jmax
ix = ix + 1
if (j .eq. jmax) goto 2
pd = cst*real(i+j-1)
c
c north
c
a(jax) = pd
if (i.eq. 1) a(jax) = a(jax)+pd
ja(jax) = ix + 1
jax = jax+1
c east
a(jax) = pd
if (j .eq. 1) a(jax) = a(jax)+pd
ja(jax) = ix + jmax
jax = jax+1
c south
2 pu = half - cst*real(i+j-3)
if ( j .gt. 1) then
a(jax) = pu
ja(jax) = ix-1
jax = jax+1
endif
c west
if ( i .gt. 1) then
a(jax) = pu
ja(jax) = ix - jmax - 1
jax = jax+1
endif
ia(ix+1) = jax
10 continue
20 continue
n = ix
return
end
|
