← Back to branch summary

~airpollution/fluidity/fluidity_airpollution

« back to all changes in this revision

Viewing changes to libalgencan/sevalus.F

  • Committer: ziyouzhj
  • Date: 2013-12-09 16:51:29 UTC
  • Revision ID: ziyouzhj@gmail.com-20131209165129-ucoetc3u0atyy05c
airpolution

Show diffs side-by-side

added added

removed removed

Lines of Context:
libalgencan/sevalus.F
 
1
C     SCALE OBJECTIVE FUNCTION AND CONSTRAINTS
 
2
 
 
3
C     ******************************************************************
 
4
C     ******************************************************************
 
5
 
 
6
      subroutine sinip(n,x,l,u,m,lambda,equatn,linear,coded,checkder,
 
7
     +inform)
 
8
 
 
9
      implicit none
 
10
 
 
11
C     SCALAR ARGUMENTS
 
12
      logical checkder
 
13
      integer inform,m,n
 
14
 
 
15
C     ARRAY ARGUMENTS
 
16
      logical coded(10),equatn(m),linear(m)
 
17
      double precision l(n),lambda(m),u(n),x(n)
 
18
 
 
19
#include "dim.par"
 
20
#include "outtyp.com"
 
21
#include "scaling.com"
 
22
#include "algparam.com"
 
23
 
 
24
C     LOCAL SCALARS
 
25
      integer i,fun,j,jcnnz,nbds,neq,file50_unit
 
26
      double precision scmax
 
27
      logical connected
 
28
 
 
29
C     LOCAL ARRAYS
 
30
      integer jcfun(jcnnzmax),jcvar(jcnnzmax)
 
31
      double precision g(nmax),jcval(jcnnzmax)
 
32
 
 
33
      neq = 0
 
34
      do j = 1,m
 
35
          if ( equatn(j) ) neq = neq + 1
 
36
      end do
 
37
 
 
38
      nbds = 0
 
39
      do i = 1,n
 
40
          if ( l(i) .gt. - 1.0d+20 ) nbds = nbds + 1
 
41
          if ( u(i) .lt.   1.0d+20 ) nbds = nbds + 1
 
42
      end do
 
43
 
 
44
      if ( iprintctl(2) ) then
 
45
          write(* ,100) n,neq,m-neq,nbds
 
46
          write(file10_unit,100) n,neq,m-neq,nbds
 
47
      end if
 
48
 
 
49
      call tinip(n,x,l,u,m,lambda,equatn,linear,coded,checkder,inform)
 
50
      if ( inform .lt. 0 ) return
 
51
 
 
52
C     Write classification line of final model
 
53
 
 
54
      if ( iprintctl(6) ) then
 
55
          file50_unit_loop: do file50_unit=10, 99
 
56
    
 
57
            inquire(unit=file50_unit, opened=connected)
 
58
    
 
59
            if (.not.connected) exit file50_unit_loop
 
60
    
 
61
          end do file50_unit_loop
 
62
          
 
63
          open(file50_unit,file='class-tabline.out')
 
64
          write(file50_unit,400) n,neq,m-neq,nbds
 
65
          close(file50_unit)
 
66
      end if
 
67
 
 
68
C     Scaling
 
69
 
 
70
      usf = 1.0d0
 
71
      do j = 1,m
 
72
          usc(j) = 1.0d0
 
73
      end do
 
74
 
 
75
      if ( scale ) then
 
76
 
 
77
          if ( m .eq. 0 ) then
 
78
              sf = 1.0d0
 
79
              if ( iprintctl(2) ) then
 
80
                  write(* ,200) 1.0d0 / sf
 
81
                  write(file10_unit,200) 1.0d0 / sf
 
82
              end if
 
83
 
 
84
              return
 
85
          end if
 
86
 
 
87
          call tsetp(n,x)
 
88
 
 
89
          if ( gjaccoded ) then
 
90
 
 
91
              call tevalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,inform)
 
92
              if ( inform .lt. 0 ) return
 
93
 
 
94
C             Scale constraints
 
95
 
 
96
              do j = 1,m
 
97
                  sc(j) = 1.0d0
 
98
              end do
 
99
 
 
100
              do i = 1,jcnnz
 
101
                  fun = jcfun(i)
 
102
 
 
103
                  sc(fun) = max( sc(fun), abs( jcval(i) ) )
 
104
              end do
 
105
 
 
106
          else
 
107
 
 
108
              call tevalg(n,x,g,inform)
 
109
              if ( inform .lt. 0 ) return
 
110
 
 
111
C             Scale constraints
 
112
 
 
113
              do j = 1,m
 
114
                  sc(j) = 1.0d0
 
115
 
 
116
                  call tevaljac(n,x,j,jcvar,jcval,jcnnz,inform)
 
117
                  if ( inform .lt. 0 ) return
 
118
 
 
119
                  do i = 1,jcnnz
 
120
                      sc(j) = max( sc(j), abs( jcval(i) ) )
 
121
                  end do
 
122
              end do
 
123
 
 
124
          end if
 
125
 
 
126
C         Scale objective function
 
127
 
 
128
          sf = 1.0d0
 
129
          do i = 1,n
 
130
              sf = max( sf, abs( g(i) ) )
 
131
          end do
 
132
 
 
133
C         Report scaling factors
 
134
 
 
135
          scmax = 0.0d0
 
136
          do j = 1,m
 
137
              scmax = max( scmax, sc(j) )
 
138
          end do
 
139
 
 
140
          if ( iprintctl(2) ) then
 
141
              write(* ,300) 1.0d0 / sf,1.0d0 / scmax
 
142
              write(file10_unit,300) 1.0d0 / sf,1.0d0 / scmax
 
143
          end if
 
144
      end if
 
145
 
 
146
C     NON-EXECUTABLE STATEMENTS
 
147
 
 
148
 100  format(/,1X,'Number of variables               : ',I7,
 
149
     +       /,1X,'Number of equality constraints    : ',I7,
 
150
     +       /,1X,'Number of inequality constraints  : ',I7,
 
151
     +       /,1X,'Number of bound constraints       : ',I7)
 
152
 
 
153
 200  format(/,1X,'Objective function scale factor   : ',1P,D7.1,
 
154
     +       /,1X,'The scaling feature was mainly developed for ',
 
155
     +            'constrained problems. For',/,1X,'unconstrained and ',
 
156
     +            'bound-constrained problem, please, set the ',
 
157
     +            'optimality',/,1X,'tolerance (related to the ',
 
158
     +            'sup-norm of the projected gradient of the',/,1X,
 
159
     +            'objective function) with a convenient value.')
 
160
 
 
161
 300  format(/,1X,'Objective function scale factor   : ',1P,D7.1,
 
162
     +       /,1X,'Smallest constraints scale factor : ',1P,D7.1)
 
163
 
 
164
 400  format(  1X,I6,1X,I6,1X,I6,1X,I6)
 
165
 
 
166
      end
 
167
 
 
168
C     ******************************************************************
 
169
C     ******************************************************************
 
170
 
 
171
      subroutine sendp(n,x,l,u,m,lambda,equatn,linear,inform)
 
172
 
 
173
      implicit none
 
174
 
 
175
C     SCALAR ARGUMENTS
 
176
      integer inform,m,n
 
177
 
 
178
C     ARRAY ARGUMENTS
 
179
      logical equatn(m),linear(m)
 
180
      double precision l(n),lambda(m),u(n),x(n)
 
181
 
 
182
#include "dim.par"
 
183
#include "scaling.com"
 
184
 
 
185
C     LOCAL SCALARS
 
186
      integer i
 
187
 
 
188
      if ( scale ) then
 
189
          do i = 1,m
 
190
              lambda(i) = lambda(i) * sf / sc(i)
 
191
          end do
 
192
 
 
193
          scale = .false.
 
194
      end if
 
195
 
 
196
      call tendp(n,x,l,u,m,lambda,equatn,linear,inform)
 
197
      if ( inform .lt. 0 ) return
 
198
 
 
199
      end
 
200
 
 
201
C     ******************************************************************
 
202
C     ******************************************************************
 
203
 
 
204
      subroutine sevalf(n,x,f,inform)
 
205
 
 
206
      implicit none
 
207
 
 
208
C     SCALAR ARGUMENTS
 
209
      integer inform,n
 
210
      double precision f
 
211
 
 
212
C     ARRAY ARGUMENTS
 
213
      double precision x(n)
 
214
 
 
215
#include "dim.par"
 
216
#include "scaling.com"
 
217
 
 
218
      call tevalf(n,x,f,inform)
 
219
      if ( inform .lt. 0 ) return
 
220
 
 
221
      if ( scale ) f = f / sf
 
222
 
 
223
      end
 
224
 
 
225
C     ******************************************************************
 
226
C     ******************************************************************
 
227
 
 
228
      subroutine sevalg(n,x,g,inform)
 
229
 
 
230
      implicit none
 
231
 
 
232
C     SCALAR ARGUMENTS
 
233
      integer inform,n
 
234
 
 
235
C     ARRAY ARGUMENTS
 
236
      double precision g(n),x(n)
 
237
 
 
238
#include "dim.par"
 
239
#include "scaling.com"
 
240
 
 
241
C     LOCAL SCALARS
 
242
      integer i
 
243
 
 
244
      call tevalg(n,x,g,inform)
 
245
      if ( inform .lt. 0 ) return
 
246
 
 
247
      if ( scale ) then
 
248
          do i = 1,n
 
249
              g(i) = g(i) / sf
 
250
          end do
 
251
      end if
 
252
 
 
253
      end
 
254
 
 
255
C     ******************************************************************
 
256
C     ******************************************************************
 
257
 
 
258
      subroutine sevalh(n,x,hlin,hcol,hval,hnnz,inform)
 
259
 
 
260
      implicit none
 
261
 
 
262
C     SCALAR ARGUMENTS
 
263
      integer inform,n,hnnz
 
264
 
 
265
C     ARRAY ARGUMENTS
 
266
      integer hcol(*),hlin(*)
 
267
      double precision hval(*),x(n)
 
268
 
 
269
#include "dim.par"
 
270
#include "scaling.com"
 
271
 
 
272
C     LOCAL SCALARS
 
273
      integer i
 
274
 
 
275
      call tevalh(n,x,hlin,hcol,hval,hnnz,inform)
 
276
      if ( inform .lt. 0 ) return
 
277
 
 
278
      if ( scale ) then
 
279
          do i = 1,hnnz
 
280
              hval(i) = hval(i) / sf
 
281
          end do
 
282
      end if
 
283
 
 
284
      end
 
285
 
 
286
C     ******************************************************************
 
287
C     ******************************************************************
 
288
 
 
289
      subroutine sevalc(n,x,ind,c,inform)
 
290
 
 
291
      implicit none
 
292
 
 
293
C     SCALAR ARGUMENTS
 
294
      integer ind,inform,n
 
295
      double precision c
 
296
 
 
297
C     ARRAY ARGUMENTS
 
298
      double precision x(n)
 
299
 
 
300
#include "dim.par"
 
301
#include "scaling.com"
 
302
 
 
303
      call tevalc(n,x,ind,c,inform)
 
304
      if ( inform .lt. 0 ) return
 
305
 
 
306
      if ( scale ) c = c / sc(ind)
 
307
 
 
308
      end
 
309
 
 
310
C     ******************************************************************
 
311
C     ******************************************************************
 
312
 
 
313
      subroutine sevaljac(n,x,ind,jcvar,jcval,jcnnz,inform)
 
314
 
 
315
      implicit none
 
316
 
 
317
C     SCALAR ARGUMENTS
 
318
      integer inform,ind,n,jcnnz
 
319
 
 
320
C     ARRAY ARGUMENTS
 
321
      integer jcvar(n)
 
322
      double precision x(n),jcval(n)
 
323
 
 
324
#include "dim.par"
 
325
#include "scaling.com"
 
326
 
 
327
C     LOCAL SCALARS
 
328
      integer i
 
329
 
 
330
      call tevaljac(n,x,ind,jcvar,jcval,jcnnz,inform)
 
331
      if ( inform .lt. 0 ) return
 
332
 
 
333
      if ( scale ) then
 
334
          do i = 1,jcnnz
 
335
              jcval(i) = jcval(i) / sc(ind)
 
336
          end do
 
337
      end if
 
338
 
 
339
      end
 
340
 
 
341
C     ******************************************************************
 
342
C     ******************************************************************
 
343
 
 
344
      subroutine sevalhc(n,x,ind,hlin,hcol,hval,hnnz,inform)
 
345
 
 
346
      implicit none
 
347
 
 
348
C     SCALAR ARGUMENTS
 
349
      integer inform,ind,n,hnnz
 
350
 
 
351
C     ARRAY ARGUMENTS
 
352
      integer hcol(*),hlin(*)
 
353
      double precision hval(*),x(n)
 
354
 
 
355
#include "dim.par"
 
356
#include "scaling.com"
 
357
 
 
358
C     LOCAL SCALARS
 
359
      integer i
 
360
 
 
361
      call tevalhc(n,x,ind,hlin,hcol,hval,hnnz,inform)
 
362
      if ( inform .lt. 0 ) return
 
363
 
 
364
      if ( scale ) then
 
365
          do i = 1,hnnz
 
366
              hval(i) = hval(i) / sc(ind)
 
367
          end do
 
368
      end if
 
369
 
 
370
      end
 
371
 
 
372
C     ******************************************************************
 
373
C     ******************************************************************
 
374
 
 
375
      subroutine sevalhl(n,x,m,lambda,hlin,hcol,hval,hnnz,inform)
 
376
 
 
377
      implicit none
 
378
 
 
379
C     SCALAR ARGUMENTS
 
380
      integer hnnz,inform,m,n
 
381
 
 
382
C     ARRAY ARGUMENTS
 
383
      integer hlin(*),hcol(*)
 
384
      double precision hval(*),lambda(m),x(n)
 
385
 
 
386
#include "dim.par"
 
387
#include "scaling.com"
 
388
 
 
389
      if ( scale ) then
 
390
          call tevalhl(n,x,m,lambda,sf,sc,hlin,hcol,hval,hnnz,inform)
 
391
          if ( inform .lt. 0 ) return
 
392
 
 
393
      else
 
394
          call tevalhl(n,x,m,lambda,usf,usc,hlin,hcol,hval,hnnz,inform)
 
395
          if ( inform .lt. 0 ) return
 
396
      end if
 
397
 
 
398
      end
 
399
 
 
400
C     ******************************************************************
 
401
C     ******************************************************************
 
402
 
 
403
      subroutine sevalhlp(n,x,m,lambda,p,hp,gothl,inform)
 
404
 
 
405
      implicit none
 
406
 
 
407
C     SCALAR ARGUMENTS
 
408
      logical gothl
 
409
      integer inform,m,n
 
410
 
 
411
C     ARRAY ARGUMENTS
 
412
      double precision hp(n),lambda(m),p(n),x(n)
 
413
 
 
414
#include "dim.par"
 
415
#include "scaling.com"
 
416
 
 
417
      if ( scale ) then
 
418
          call tevalhlp(n,x,m,lambda,sf,sc,p,hp,gothl,inform)
 
419
          if ( inform .lt. 0 ) return
 
420
 
 
421
      else
 
422
          call tevalhlp(n,x,m,lambda,usf,usc,p,hp,gothl,inform)
 
423
          if ( inform .lt. 0 ) return
 
424
      end if
 
425
 
 
426
      end
 
427
 
 
428
C     ******************************************************************
 
429
C     ******************************************************************
 
430
 
 
431
      subroutine sevalfc(n,x,f,m,c,inform)
 
432
 
 
433
      implicit none
 
434
 
 
435
C     SCALAR ARGUMENTS
 
436
      integer inform,m,n
 
437
      double precision f
 
438
 
 
439
C     ARRAY ARGUMENTS
 
440
      double precision c(m),x(n)
 
441
 
 
442
#include "dim.par"
 
443
#include "scaling.com"
 
444
 
 
445
C     LOCAL SCALARS
 
446
      integer j
 
447
 
 
448
      call tevalfc(n,x,f,m,c,inform)
 
449
      if ( inform .lt. 0 ) return
 
450
 
 
451
      if ( scale ) then
 
452
          f = f / sf
 
453
 
 
454
          do j = 1,m
 
455
              c(j) = c(j) / sc(j)
 
456
          end do
 
457
      end if
 
458
 
 
459
      end
 
460
 
 
461
C     ******************************************************************
 
462
C     ******************************************************************
 
463
 
 
464
      subroutine sevalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,inform)
 
465
 
 
466
      implicit none
 
467
 
 
468
C     SCALAR ARGUMENTS
 
469
      integer inform,jcnnz,m,n
 
470
 
 
471
C     ARRAY ARGUMENTS
 
472
      integer jcfun(*),jcvar(*)
 
473
      double precision g(n),jcval(*),x(n)
 
474
 
 
475
#include "dim.par"
 
476
#include "scaling.com"
 
477
 
 
478
C     LOCAL SCALARS
 
479
      integer i
 
480
 
 
481
      call tevalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,inform)
 
482
      if ( inform .lt. 0 ) return
 
483
 
 
484
      if ( scale ) then
 
485
          do i = 1,n
 
486
              g(i) = g(i) / sf
 
487
          end do
 
488
 
 
489
          do i = 1,jcnnz
 
490
              jcval(i) = jcval(i) / sc(jcfun(i))
 
491
          end do
 
492
      end if
 
493
 
 
494
      end
 
495
 
 
496
C     ******************************************************************
 
497
C     ******************************************************************
 
498
 
 
499
      subroutine ssetp(n,x)
 
500
 
 
501
      implicit none
 
502
 
 
503
C     SCALAR ARGUMENTS
 
504
      integer n
 
505
 
 
506
C     ARRAY ARGUMENTS
 
507
      double precision x(n)
 
508
 
 
509
      call tsetp(n,x)
 
510
 
 
511
      end
 
512
 
 
513
C     ******************************************************************
 
514
C     ******************************************************************
 
515
 
 
516
      subroutine sunsetp()
 
517
 
 
518
      implicit none
 
519
 
 
520
      call tunsetp()
 
521
 
 
522
      end