File: srqfn.f

package info (click to toggle)
r-cran-quantreg 5.38-1
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 3,080 kB
  • sloc: fortran: 7,125; makefile: 2
file content (343 lines) | stat: -rw-r--r-- 13,418 bytes parent folder | download | duplicates (5)
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
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
      subroutine srqfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
     &                  dsub,jdsub,nnzemax,e,je,ie,nsubmax,lindx,xlindx,
     &                  nnzlmax,lnz,xlnz,iw,iwmax,iwork,xsuper,tmpmax,
     &                  tmpvec,wwm,wwn,cachsz,level,x,s,u,c,y,b,small,
     &                  ierr,maxit,timewd)
      integer nnza,m,n,nnzdmax,nnzemax,iwmax,
     &        nnzlmax,nsubmax,cachsz,level,tmpmax,ierr,maxit,
     &        ja(nnza),jao(nnza),jdsub(nnzemax+1),jd(nnzdmax),
     &        ia(n+1),iao(m+1),id(m+1),lindx(nsubmax),xlindx(m+1),
     &        iw(m,5),xlnz(m+1),iwork(iwmax),xsuper(m+1),je(nnzemax),
     &        ie(m+1)
      double precision small,
     &                 a(nnza),ao(nnza),dsub(nnzemax+1),d(nnzdmax),
     &                 lnz(nnzlmax),c(n),y(m),wwm(m,3),tmpvec(tmpmax),
     &                 wwn(n,14),x(n),s(n),u(n),e(nnzemax),b(m)
      double precision timewd(7)
      call slpfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
     &          dsub,jdsub,nsubmax,lindx,xlindx,nnzlmax,lnz,
     &          xlnz,iw(1,1),iw(1,2),iwmax,iwork,iw(1,3),iw(1,4),
     &          xsuper,iw(1,5),tmpmax,tmpvec,wwm(1,2),cachsz,
     &          level,x,s,u,c,y,b,wwn(1,1),wwn(1,2),wwn(1,3),
     &          wwn(1,4),nnzemax,e,je,ie,wwm(1,3),wwn(1,5),wwn(1,6),
     &          wwn(1,7),wwn(1,8),wwn(1,9),wwn(1,10),wwn(1,11),
     &          wwn(1,12),wwn(1,13),wwn(1,14),wwm(1,1),small,ierr,
     &          maxit,timewd)
      return
      end
      subroutine slpfn(n,m,nnza,a,ja,ia,ao,jao,iao,nnzdmax,d,jd,id,
     &                dsub,jdsub,nsubmax,lindx,xlindx,nnzlmax,lnz,
     &                xlnz,invp,perm,iwmax,iwork,colcnt,snode,xsuper,
     &                split,tmpmax,tmpvec,newrhs,cachsz,level,x,s,u,
     &                c,y,b,r,z,w,q,nnzemax,e,je,ie,dy,dx,ds,dz,dw,dxdz,
     &                dsdw,xi,xinv,sinv,ww1,ww2,small,ierr,maxit,timewd)
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c Sparse implentation of LMS's interior point method via 
c    Ng-Peyton's sparse Cholesky factorization for sparse 
c    symmetric positive definite
c INPUT:
c     n -- the number of row in the coefficient matrix A'
c     m -- the number of column in the coefficient matrix A'
c     nnza -- the number of non-zero elements in A'
c     a -- an nnza-vector of non-zero values of the design 
c          matrix (A') stored in csr format 
c     ja -- an nnza-vector of indices of the non-zero elements of
c           the coefficient matrix 
c     ia -- an (n+1)-vector of pointers to the begining of each
c           row in a and ja
c     ao -- an nnza-vector of work space for the transpose of
c           the design matrix stored in csr format or the 
c           design matrix stored in csc format
c     jao -- an nnza-vector of work space for the indices of the 
c            transpose of the design matrix 
c     iao -- an (n+1)-vector of pointers to the begining of each
c            column in ao and jao
c     nnzdmax -- upper bound of the non-zero elements in AA'
c     d -- an nnzdmax-vector of non-zero values of AQ^(-1)
c     jd -- an nnzdmax-vector of indices in d
c     id -- an (m+1)-vector of pointers to the begining of each
c           row in d and jd
c     dsub -- the values of e excluding the diagonal elements
c     jdsub -- the indices to dsub
c     nsubmax -- upper bound of the dimension of lindx
c     lindx -- an nsub-vector of interger which contains, in 
c           column major order, the row subscripts of the nonzero
c           entries in L in a compressed storage format
c     xlindx -- an (m+1)-vector of integer of pointers for lindx
c     nnzlmax -- the upper bound of the non-zero entries in 
c                L stored in lnz, including the diagonal entries
c     lnz -- First contains the non-zero entries of d; later
c            contains the entries of the Cholesky factor
c     xlnz -- column pointer for L stored in lnz
c     invp -- an n-vector of integer of inverse permutation 
c             vector
c     perm -- an n-vector of integer of permutation vector
c     iw -- integer work array of length m
c     iwmax -- upper bound of the general purpose integer 
c              working storage iwork; set at 7*m+3
c     iwork -- an iwsiz-vector of integer as work space
c     colcnt -- array of length m, containing the number of 
c               non-zeros in each column of the factor, including
c               the diagonal entries
c     snode -- array of length m for recording supernode 
c              membership
c     xsuper -- array of length m+1 containing the supernode
c               partitioning
c     split -- an m-vector with splitting of supernodes so that 
c              they fit into cache
c     tmpmax -- upper bound of the dimension of tmpvec
c     tmpvec -- a tmpmax-vector of temporary vector
c     newrhs -- extra work vector for right-hand side and 
c               solution
c     cachsz -- size of the cache (in kilobytes) on the target 
c               machine
c     level -- level of loop unrolling while performing numerical
c              factorization
c     x -- an n-vector, the initial feasible solution in the primal
c              that corresponds to the design matrix A'
c     s -- an n-vector 
c     u -- an n-vector of upper bound for x
c     c -- an n-vector, usually the "negative" of
c          the pseudo response
c     y -- an m-vector, the initial dual solution 
c     b -- an n-vector, the rhs of the equality constraint
c          usually X'a = (1-tau)X'e in the rq setting
c     r -- an n-vector of residuals
c     z -- an n-vector of the dual slack variable
c     w -- an n-vector 
c     q -- an n-vector of work array containing the diagonal 
c          elements of the Q^(-1) matrix
c     nnzemax -- upper bound of the non-zero elements in AA'
c     e -- an nnzdmax-vector containing the non-zero entries of
c          AQ^(-1)A' stored in csr format
c     je -- an nnzemax-vector of indices for e
c     ie -- an (m+1)-vector of pointers to the begining of each
c           row in e and je
c     dy -- work array
c     dx -- work array
c     ds -- work array
c     dz -- work array
c     dw -- work array
c     dxdz -- work array
c     dsdw -- work arry
c     xi -- work array
c     xinv -- work array
c     sinv -- work array
c     ww1 -- work array
c     ww2 -- work array
c     small -- convergence tolerance for inetrior algorithm
c     ierr -- error flag
c       1 -- insufficient work space in call to extract
c       2 -- nnzd > nnzdmax
c       3 -- insufficient storage in iwork when calling ordmmd;
c       4 -- insufficient storage in iwork when calling sfinit;
c       5 -- nnzl > nnzlmax when calling sfinit
c       6 -- nsub > nsubmax when calling sfinit
c       7 -- insufficient work space in iwork when calling symfct
c       8 -- inconsistancy in input when calling symfct
c       9 -- tmpsiz > tmpmax when calling bfinit; increase tmpmax
c       10 -- nonpositive diagonal encountered when calling 
c            blkfct, the matrix is not positive definite
c       11 -- insufficient work storage in tmpvec when calling 
c            blkfct
c       12 -- insufficient work storage in iwork when calling 
c            blkfct
c     maxit -- upper limit of the iteration; on return holds the
c              number of iterations
c     timew -- amount of time to execute this subroutine
c OUTPUT:
c     y -- an m-vector of primal solution
c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
      integer nnza,m,n,nsuper,nnzdmax,nnzemax,iwmax,nnzd,
     &        nnzlmax,nsubmax,cachsz,level,tmpmax,ierr,maxit,it,
     &        ja(nnza),jao(nnza),jdsub(nnzemax+1),jd(nnzdmax),
     &        ia(n+1),iao(m+1),id(m+1),lindx(nsubmax),xlindx(m+1),
     &        invp(m),perm(m),xlnz(m+1),iwork(iwmax),
     &        colcnt(m),snode(m),xsuper(m+1),split(m),je(nnzemax),
     &        ie(m+1)
      double precision ddot,gap,zero,one,beta,small,deltap,deltad,mu,g,
     &                 a(nnza),ao(nnza),dsub(nnzemax+1),d(nnzdmax),
     &                 lnz(nnzlmax),c(n),b(m),newrhs(m),y(m),
     &                 tmpvec(tmpmax),r(n),z(n),w(n),x(n),s(n),
     &                 u(n),q(n),
     &                 e(nnzemax),dy(m),dx(n),ds(n),dz(n),dw(n),
     &                 dxdz(n),dsdw(n),xinv(n),sinv(n),xi(n),
     &                 ww1(n),ww2(m)
      double precision timewd(7)
      real gtimer,timbeg,timend
      external smxpy1,smxpy2,smxpy4,smxpy8
      external mmpy1,mmpy2,mmpy4,mmpy8
      parameter (beta=9.995d-1, one=1.0d0, zero=0.0d0)
      do i = 1,7
         timewd(i) = 0.0
      enddo
      it = 0
      nnzd = ie(m+1) - 1
      nnzdsub = nnzd - m
c
c  Compute the initial gap
c
       gap = ddot(n,z,1,x,1) + ddot(n,w,1,s,1)
c
c  Start iteration
c
   20 continue
      if(gap .lt. small .or. it .gt. maxit) goto 30
      it = it + 1
c
c  Create the diagonal matrix Q^(-1) stored in q as an n-vector
c  and update the residuals in r
c
      do i=1,n
         q(i) = one/(z(i)/x(i)+w(i)/s(i))
         r(i) = z(i) - w(i)
      enddo
c
c  Obtain AQ^(-1) and store in d,jd,id in csr format
c
      call amudia(m,1,ao,jao,iao,q,d,jd,id)
c
c  Obtain AQ^(-1)A' and store in e,je,ie in csr format
c
      call amub(m,m,1,d,jd,id,a,ja,ia,e,je,ie,nnzemax,iwork,ierr)
      if (ierr .ne. 0) then
         ierr = 2
         go to 100
      endif
c
c  Extract the non-diagonal structure of e,je,ie and store in dsub,jdsub
c
      call extract(e,je,ie,dsub,jdsub,m,nnzemax,nnzemax+1,ierr)
      if (ierr .ne. 0) then
         ierr = 1
         go to 100
      endif
c
c Compute b - Ax + AQ^(-1)r and store it in c in two steps
c  First: store Ax in ww2
      call amux(m,x,ww2,ao,jao,iao)
c
c  Second: save AQ^(-1)r in c temporarily
c
      call amux(m,r,c,d,jd,id)
      do i = 1,m
         c(i) = b(i) - ww2(i) + c(i)
      enddo
c
c  Compute dy = (AQ^(-1)A')^(-1)(b-Ax+AQ^(-1)r); result returned via dy
c
c Call chlfct to perform Cholesky's decomposition of e,je,ie
c
      call chlfct(m,xlindx,lindx,invp,perm,iwork,nnzdsub,jdsub,
     &            colcnt,nsuper,snode,xsuper,nnzlmax,nsubmax,xlnz,lnz,
     &            ie,je,e,cachsz,tmpmax,level,tmpvec,split,ierr,it,
     &            timewd)
      if (ierr .ne. 0) go to 100
c
c Call blkslv: Numerical solution for the new rhs stored in b
c
      do i = 1,m
         newrhs(i) = c(perm(i))
      enddo
      timbeg = gtimer()
      call blkslv(nsuper,xsuper,xlindx,lindx,xlnz,lnz,newrhs)
      timend = gtimer()
      timewd(7) = timewd(7) + timend - timbeg
      do i = 1,m
         dy(i) = newrhs(invp(i))
      enddo
c
c  Compute dx = Q^(-1)(A'dy - r), ds = -dx, dz  and dw
c
      call amux(n,dy,dx,a,ja,ia)
      do i=1,n
         dx(i) = q(i) * (dx(i) - r(i))
         ds(i) = -dx(i)
         dz(i) = -z(i) * (one + dx(i) / x(i))
         dw(i) = -w(i) * (one + ds(i) / s(i))
      enddo
c
c Compute the maximum allowable step lengths
c
      call bound(x,dx,s,ds,z,dz,w,dw,n,beta,deltap,deltad)
      if (deltap * deltad .lt. one) then
c
c Update mu
c
         mu = ddot(n,z,1,x,1) + ddot(n,w,1,s,1)
         g = ddot(n,z,1,x,1) + deltap*ddot(n,z,1,dx,1)
     &       + deltad*ddot(n,dz,1,x,1) + deltad*deltap*ddot(n,dz,1,dx,1)
     &       + ddot(n,w,1,s,1) + deltap*ddot(n,w,1,ds,1)
     &       + deltad*ddot(n,dw,1,s,1) + deltad*deltap*ddot(n,dw,1,ds,1)
         mu = mu*((g/mu)**3)/(2.d0*dfloat(n))
c
c Compute dxdz and dsdw
c
         do i = 1,n
             dxdz(i) = dx(i)*dz(i)
             dsdw(i) = ds(i)*dw(i)
             xinv(i) = one/x(i)
             sinv(i) = one/s(i)
             xi(i) = xinv(i) * dxdz(i) - sinv(i) * dsdw(i)
     &               - mu * (xinv(i) - sinv(i))
             ww1(i) = q(i) * xi(i)
         enddo
c
c Compute AQ^(-1)(dxdz - dsdw - mu(X^(-1) - S^(-1))) and
c store it in ww2 temporarily
c
         call amux(m,ww1,ww2,ao,jao,iao)
         do i = 1,m
            c(i) = c(i) + ww2(i)
         enddo
c
c
c Compute dy and return the result in dy
c
c Call blkslv: Numerical solution for the new rhs stored in b
c
      do i = 1,m
         newrhs(i) = c(perm(i))
      enddo
      timbeg = gtimer()
      call blkslv(nsuper,xsuper,xlindx,lindx,xlnz,lnz,newrhs)
      timend = gtimer()
      timewd(7) = timewd(7) + timend - timbeg
      do i = 1,m
         dy(i) = newrhs(invp(i))
      enddo
c
c  Compute dx = Q^(-1)(A'dy - r + mu(X^(-1) - S^(-1)) -dxdz + dsdw),
c  ds = -dx, dz  and dw
c
         call amux(n,dy,dx,a,ja,ia)
         do i=1,n
            dx(i) = q(i) * (dx(i) - xi(i) - r(i))
            ds(i) = -dx(i)
            dz(i) = -z(i) + xinv(i)*(mu - z(i)*dx(i) - dxdz(i))
            dw(i) = -w(i) + sinv(i)*(mu - w(i)*ds(i) - dsdw(i))
         enddo
c
c Compute the maximum allowable step lengths
c
         call bound(x,dx,s,ds,z,dz,w,dw,n,beta,deltap,deltad)
      endif
c
c Take the step
c
      call daxpy(n,deltap,dx,1,x,1)
      call daxpy(n,deltap,ds,1,s,1)
      call daxpy(n,deltad,dw,1,w,1)
      call daxpy(n,deltad,dz,1,z,1)
      call daxpy(m,deltad,dy,1,y,1)
      gap = ddot(n,z,1,x,1) + ddot(n,w,1,s,1)
      goto 20
   30 continue
  100 continue
      maxit = it
      return
      end