← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/sparse/lspops.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/sparse/lspops.f
 
1
      subroutine lspops
 
2
c     
 
3
c     operations on boolean sparse matrices
 
4
c     
 
5
c     Copyright INRIA
 
6
      include '../stack.h'
 
7
      integer op
 
8
c     
 
9
      integer iadr,sadr
 
10
c     
 
11
      integer star,dstar,dot,colon
 
12
      integer less,great,equal,ou,et,non
 
13
      integer insert,extrac
 
14
      integer top0
 
15
      logical isany
 
16
c
 
17
      data star/47/,dstar/62/,dot/51/,colon/44/
 
18
      data less/59/,great/60/,equal/50/
 
19
      data ou/57/,et/58/,non/61/
 
20
      data insert/2/,extrac/3/
 
21
c     
 
22
      iadr(l)=l+l-1
 
23
      sadr(l)=(l/2)+1
 
24
c     
 
25
      op=fin
 
26
c     
 
27
      if (ddt .eq. 4) then
 
28
         write(buf(1:4),'(i4)') fin
 
29
         call basout(io,wte,' lspops op: '//buf(1:4))
 
30
      endif
 
31
c     
 
32
      top0=top
 
33
      lw=lstk(top+1)+1
 
34
      if(op.eq.extrac) goto 70
 
35
      if(op.eq.insert) goto 80
 
36
 
 
37
      it2=0
 
38
      goto (04,03,02,01) rhs
 
39
      call error(39)
 
40
      return
 
41
c     
 
42
 01   il4=iadr(lstk(top))
 
43
      if(istk(il4).lt.0) il4=iadr(istk(il4+1))
 
44
      m4=istk(il4+1)
 
45
      n4=istk(il4+2)
 
46
      it4=istk(il4+3)
 
47
      if(istk(il4).eq.6) then
 
48
         nel4=istk(il4+4)
 
49
         irc4=il4+5
 
50
         l4=sadr(irc4+m4+nel4)
 
51
      else
 
52
         nel4=m4*n4
 
53
         l4=sadr(il4+4)
 
54
      endif
 
55
      mn4=m4*n4
 
56
      top=top-1
 
57
c     
 
58
 02   il3=iadr(lstk(top))
 
59
      if(istk(il3).lt.0) il3=iadr(istk(il3+1))
 
60
      m3=istk(il3+1)
 
61
      n3=istk(il3+2)
 
62
      it3=istk(il3+3)
 
63
      if(istk(il3).eq.6) then
 
64
         nel3=istk(il3+4)
 
65
         irc3=il3+5
 
66
         l3=sadr(irc3+m3+nel3)
 
67
      else
 
68
         l3=sadr(il3+4)
 
69
         nel3=m3*n3
 
70
      endif
 
71
      mn3=m3*n3
 
72
      top=top-1
 
73
c     
 
74
 03   il2=iadr(lstk(top))
 
75
      if(istk(il2).lt.0) il2=iadr(istk(il2+1))
 
76
      m2=istk(il2+1)
 
77
      n2=istk(il2+2)
 
78
      it2=istk(il2+3)
 
79
      if(istk(il2).eq.6) then
 
80
         nel2=istk(il2+4)
 
81
         irc2=il2+5
 
82
         l2=sadr(irc2+m2+nel2)
 
83
      else
 
84
         l2=sadr(il2+4)
 
85
         nel2=m2*n2
 
86
      endif
 
87
      mn2=m2*n2
 
88
      top=top-1
 
89
c     
 
90
 04   il1=iadr(lstk(top))
 
91
      il1r=il1
 
92
      if(istk(il1).lt.0) il1=iadr(istk(il1+1))
 
93
      m1=istk(il1+1)
 
94
      n1=istk(il1+2)
 
95
      it1=istk(il1+3)
 
96
      if(istk(il1).eq.6) then
 
97
         nel1=istk(il1+4)
 
98
         irc1=il1+5
 
99
         l1=sadr(irc1+m1+nel1)
 
100
      else
 
101
         l1=sadr(il1+4)
 
102
         nel1=m1*n1
 
103
      endif
 
104
      mn1=m1*n1
 
105
      top=top-1
 
106
c     
 
107
c     operations binaires et ternaires
 
108
c     --------------------------------
 
109
c     
 
110
      top=top+1
 
111
      itr=max(it1,it2)
 
112
c     
 
113
      fun = 0
 
114
c     
 
115
c     column concatenation
 
116
      if(op.eq.1) goto 65
 
117
c     row concatenation
 
118
      if(op.eq.4) goto 66
 
119
c           :  +  -  * /  \  =          '
 
120
      goto(07,07,07,07,07,07,130,05,05,60) op+1-colon
 
121
      if(op.eq.ou.or.op.eq.et) goto 20
 
122
      if(op.eq.non) goto 30
 
123
c     
 
124
 05   if(op.eq.dstar) goto 07
 
125
      if(op.ge.3*dot+star) goto 07
 
126
      if(op.ge.2*dot+star) goto 07
 
127
      if(op.ge.less+equal) goto 130
 
128
      if(op.ge.dot+star) goto 07
 
129
      if(op.ge.less) goto 130
 
130
 
 
131
 06   call error(43)
 
132
      return
 
133
c
 
134
 07   fin=-fin
 
135
      top=top0
 
136
      go to 999
 
137
c     
 
138
c     ou/et logique
 
139
 20   if(istk(il1).ne.6.or.istk(il2).ne.6) then
 
140
         fin=-fin
 
141
         top=top0
 
142
         go to 999
 
143
      endif
 
144
      if(mn1.eq.1.and.mn2.gt.1) then
 
145
         top=top0
 
146
         fin=-fin
 
147
         return
 
148
      elseif(mn2.eq.1.and.mn1.gt.1) then
 
149
         top=top0
 
150
         fin=-fin
 
151
         return
 
152
      else if (n1 .ne. n2.or.m1.ne.m2) then
 
153
         call error(60)
 
154
         return
 
155
      endif
 
156
      irc=iadr(lw)
 
157
      nelmx=iadr(lstk(bot))-irc-m1-10
 
158
      lw=sadr(irc+m1+nelmx)
 
159
      err=lw-lstk(bot)
 
160
      if(err.gt.0) then
 
161
         call error(17)
 
162
         return
 
163
      endif  
 
164
      nel=nelmx
 
165
      if(fin.eq.ou) then
 
166
         call lspasp(m1,n1,nel1,istk(irc1),nel2,istk(irc2),nel,
 
167
     $        istk(irc),ierr)
 
168
      else
 
169
         call lspxsp(m1,n1,nel1,istk(irc1),nel2,istk(irc2),nel,
 
170
     $        istk(irc),ierr)
 
171
      endif
 
172
      if(ierr.ne.0) then
 
173
         call error(17)
 
174
         return
 
175
      endif
 
176
      istk(il1+3)=0
 
177
      istk(il1+4)=nel
 
178
      call icopy(m1+nel,istk(irc),1,istk(irc1),1)
 
179
      l1=sadr(irc1+m1+nel)
 
180
      lstk(top+1)=l1
 
181
      go to 999
 
182
 
 
183
c NOT
 
184
 30   continue
 
185
      lw=iadr(lstk(top+1))
 
186
      err=sadr(lw+m1*n1-nel1+m1)-lstk(bot)
 
187
      if(err.gt.0) then
 
188
         call error(17)
 
189
         return
 
190
      endif
 
191
 
 
192
      istk(il1+4)=m1*n1-nel1
 
193
      ij2=lw
 
194
      ij1=irc1+m1
 
195
c     may be improved
 
196
      do 35 i=0,m1-1
 
197
         do 31 j=1,n1
 
198
            istk(ij2+j-1)=j
 
199
 31      continue
 
200
         do 32 j=1,istk(irc1+i)
 
201
            istk(ij2+istk(ij1+j-1)-1)=0
 
202
 32      continue
 
203
         j1=0
 
204
         do 33 j=1,n1
 
205
            if(istk(ij2+j-1).ne.0) then
 
206
               j1=j1+1
 
207
               istk(ij2+j1-1)=istk(ij2+j-1)
 
208
            endif
 
209
 33      continue
 
210
         ij2=ij2+j1
 
211
         ij1=ij1+istk(irc1+i)
 
212
         istk(irc1+i)=n1-istk(irc1+i)
 
213
 35   continue
 
214
      call icopy(ij2-lw,istk(lw),1,istk(irc1+m1),1)
 
215
      lstk(top+1)=sadr(irc1+m1+ij2-lw)
 
216
      goto 999
 
217
 
 
218
c     
 
219
c     
 
220
c     transposition
 
221
 60   istk(il1+1)=n1
 
222
      istk(il1+2)=m1
 
223
      if(nel1.eq.0) then
 
224
         lw=sadr(il1+5+n1)
 
225
         err=lw-lstk(bot)
 
226
         if(err.gt.0) then
 
227
            call error(17)
 
228
            return
 
229
         endif
 
230
         call iset(n1,0,istk(il1+5),1)
 
231
         lstk(top+1)=lw
 
232
         goto 999
 
233
      endif
 
234
      ia=iadr(lw)
 
235
      iat=ia+m1+1
 
236
      irc=iat+n1+1
 
237
      lw=sadr(irc+n1+nel1)
 
238
      err=lw-lstk(bot)
 
239
      if(err.gt.0) then
 
240
         call error(17)
 
241
         return
 
242
      endif
 
243
      istk(ia)=1
 
244
      do 61 i=1,m1
 
245
         istk(ia+i)=istk(ia+i-1)+istk(irc1+i-1)
 
246
 61   continue
 
247
      call lspt(m1,n1,nel1,istk(irc1),istk(ia),
 
248
     $     istk(iat),istk(irc))
 
249
      call icopy(n1+nel1,istk(irc),1,istk(irc1),1)
 
250
      l1=sadr(irc1+n1+nel1)
 
251
      lstk(top+1)=l1
 
252
      goto 999
 
253
c     
 
254
c     concatenation [a b]
 
255
 65   continue
 
256
      if(m1.lt.0.or.m2.lt.0) then
 
257
         call error(14)
 
258
         return
 
259
      endif
 
260
      if(m2.eq.0) then
 
261
         return
 
262
      elseif(m1.eq.0) then
 
263
         call icopy(5+m2+nel2,istk(il2),1,istk(il1),1)
 
264
         l1=sadr(il1+5+m2+nel2)
 
265
         lstk(top+1)=l1
 
266
         return
 
267
      elseif(m1.ne.m2) then
 
268
         call error(5)
 
269
         return
 
270
      endif
 
271
      if(istk(il1).ne.6.or.istk(il2).ne.6) then
 
272
         top=top0
 
273
         fin=-fin
 
274
         return
 
275
      endif
 
276
c
 
277
      nelr=nel1+nel2
 
278
      istk(il1+2)=n1+n2
 
279
      istk(il1+3)=0
 
280
      istk(il1+4)=nelr
 
281
      irc=iadr(lw)
 
282
      lw=sadr(irc+m1+nelr)
 
283
      err=lw-lstk(bot)
 
284
      if(err.gt.0) then
 
285
         call error(17)
 
286
         return
 
287
      endif
 
288
      call lspcsp(0,m1,n1,nel1,istk(irc1),
 
289
     $     m2,n2,nel2,istk(irc2),
 
290
     $     nelr,istk(irc))
 
291
      call icopy(m1+nelr,istk(irc),1,istk(irc1),1)
 
292
      l1=sadr(irc1+m1+nelr)
 
293
      lstk(top+1)=l1
 
294
      return
 
295
c     
 
296
c     concatenation [a;b]
 
297
 66   continue
 
298
      if(n1.lt.0.or.n2.lt.0) then
 
299
         call error(14)
 
300
         return
 
301
      endif
 
302
      if(n2.eq.0) then
 
303
         goto 999
 
304
      elseif(n1.eq.0)then
 
305
         call icopy(5+m2+nel2,istk(il2),1,istk(il1),1)
 
306
         l1=sadr(il1+5+m2+nel2)
 
307
         lstk(top+1)=l1
 
308
         goto 999
 
309
      elseif(n1.ne.n2) then
 
310
         call error(6)
 
311
         return
 
312
      endif
 
313
      if(istk(il1).ne.6.or.istk(il2).ne.6) then
 
314
         top=top0
 
315
         fin=-fin
 
316
         return
 
317
      endif
 
318
 
 
319
 
 
320
      nelr=nel1+nel2
 
321
      istk(il1+1)=m1+m2
 
322
      istk(il1+3)=0
 
323
      istk(il1+4)=nelr
 
324
      irc=iadr(lw)
 
325
      lw=sadr(irc+m1+m2+nelr)
 
326
      err=lw-lstk(bot)
 
327
      if(err.gt.0) then
 
328
         call error(17)
 
329
         return
 
330
      endif
 
331
      call lspcsp(1,m1,n1,nel1,istk(irc1),
 
332
     $        m2,n2,nel2,istk(irc2),
 
333
     $        nelr,istk(irc))
 
334
      call icopy(m1+m2+nelr,istk(irc),1,istk(irc1),1)
 
335
      l1=sadr(irc1+m1+m2+nelr)
 
336
      lstk(top+1)=l1
 
337
      goto 999
 
338
c     
 
339
c     extraction
 
340
c     
 
341
 70   continue
 
342
      if(rhs.gt.2) goto 75
 
343
c     arg2(arg1)
 
344
c     get arg2
 
345
      il2=iadr(lstk(top))
 
346
      if(istk(il2).lt.0) il2=iadr(istk(il2+1))
 
347
      m2=istk(il2+1)
 
348
      n2=istk(il2+2)
 
349
      it2=istk(il2+3)
 
350
      nel2=istk(il2+4)
 
351
      irc2=il2+5
 
352
      l2=sadr(irc2+m2+nel2)
 
353
      mn2=m2*n2
 
354
      top=top-1
 
355
c     get arg1
 
356
      il1=iadr(lstk(top))
 
357
      ilrs=il1
 
358
      if(istk(il1).lt.0) il1=iadr(istk(il1+1))
 
359
      m1=istk(il1+1)
 
360
c
 
361
      if(mn2.eq.0) then 
 
362
c     .  arg2=[]
 
363
         ilrs=iadr(lstk(top))
 
364
         istk(ilrs)=1
 
365
         istk(ilrs+1)=0
 
366
         istk(ilrs+2)=0
 
367
         istk(ilrs+3)=0
 
368
         lstk(top+1)=sadr(ilrs+4)+1
 
369
         goto 999
 
370
      elseif(m2.lt.0) then
 
371
c     .  arg2=eye
 
372
         call error(14)
 
373
         return
 
374
      elseif(m1.lt.0) then
 
375
c     .  arg2(:), just reshape to column vector
 
376
         if(n2.eq.1) then
 
377
            call icopy(5+m2+nel2,istk(il2),1,istk(ilrs),1)
 
378
            l1=sadr(ilrs+5+m2+nel2)
 
379
            lstk(top+1)=l1
 
380
         else
 
381
c     .     reshape to column vector
 
382
            ilrs=iadr(lstk(top))
 
383
            istk(ilrs)=6
 
384
            istk(ilrs+1)=mn2
 
385
            istk(ilrs+2)=1
 
386
            istk(ilrs+3)=0
 
387
            istk(ilrs+4)=nel2
 
388
            irc1=ilrs+5
 
389
            l1=sadr(ilrs+5+m2*n2+nel2)
 
390
 
 
391
            ircr=iadr(lw)
 
392
            iw=ircr+m2*n2+nel2
 
393
            lw=sadr(iw+3*nel2)
 
394
            err=lw-lstk(bot)
 
395
            if(err.gt.0) then
 
396
               call error(17)
 
397
               return 
 
398
            endif
 
399
            call lspmat(m2,n2,nel2,istk(irc2),m2*n2,istk(ircr),istk(iw))
 
400
            call icopy(m2*n2+nel2,istk(ircr),1,istk(irc1),1)
 
401
            lstk(top+1)=l1
 
402
         endif
 
403
         return
 
404
      elseif(m2.gt.1.and.n2.gt.1) then
 
405
c     .  call macro coded operation
 
406
         top=top0
 
407
         fin=-fin
 
408
         return
 
409
      endif
 
410
c     check and convert indices variable
 
411
      call indxg(il1,mn2,ilr,mi,mx,lw,1)
 
412
      if(err.gt.0) return
 
413
      if(mx.gt.mn2) then
 
414
         call error(21)
 
415
         return
 
416
      endif
 
417
 72   if(mi.eq.0) then
 
418
c     arg2([])
 
419
         ilrs=iadr(lstk(top))
 
420
         istk(ilrs)=1
 
421
         istk(ilrs+1)=0
 
422
         istk(ilrs+2)=0
 
423
         istk(ilrs+3)=0
 
424
         lstk(top+1)=sadr(ilrs+4)+1
 
425
         goto 999
 
426
      endif
 
427
c     set output sizes
 
428
      if (m2 .gt. 1.or.m1.lt.0) then
 
429
c     .  column vector
 
430
         m=mi
 
431
         n=-1
 
432
         mr = mi
 
433
         nr = 1
 
434
      else
 
435
c     .  row vector
 
436
         m=-1
 
437
         n=mi
 
438
         nr = mi
 
439
         mr = 1
 
440
      endif
 
441
      lptr=iadr(lw)
 
442
      irc=lptr+m2+1
 
443
      lw=sadr(irc+mr)
 
444
      nelr=iadr(lstk(bot))-iadr(lw)
 
445
      if(nelr.le.0) then
 
446
         err=lw-lstk(bot)
 
447
         call error(17)
 
448
         return
 
449
      endif
 
450
      lw=sadr(irc+mr+nelr)
 
451
      nel=nelr
 
452
 
 
453
      call lspe2(m2,n2,nel2,istk(irc2),
 
454
     $        istk(ilr),m,istk(ilr),n,mr,nr,
 
455
     $        nelr,istk(irc),istk(lptr),ierr)
 
456
      ilrs=iadr(lstk(top))
 
457
      istk(ilrs)=6
 
458
      istk(ilrs+1)=mr
 
459
      istk(ilrs+2)=nr
 
460
      istk(ilrs+3)=0
 
461
      istk(ilrs+4)=nelr
 
462
      call icopy(m+nelr,istk(irc),1,istk(ilrs+5),1)
 
463
      l1=sadr(ilrs+5+mr+nelr)
 
464
      lstk(top+1)=l1
 
465
      go to 999
 
466
 
 
467
c     
 
468
c     arg3(arg1,arg2)
 
469
 75   if(rhs.gt.3) then
 
470
         call error(36)
 
471
         return
 
472
      endif
 
473
c     get arg3
 
474
      il3=iadr(lstk(top))
 
475
      if(istk(il3).lt.0) il3=iadr(istk(il3+1))
 
476
      m3=istk(il3+1)
 
477
      n3=istk(il3+2)
 
478
      it3=istk(il3+3)
 
479
      nel3=istk(il3+4)
 
480
      irc3=il3+5
 
481
      l3=sadr(irc3+m3+nel3)
 
482
      mn3=m3*n3
 
483
      top=top-1
 
484
c     get arg2
 
485
      il2=iadr(lstk(top))
 
486
      if(istk(il2).lt.0) il2=iadr(istk(il2+1))
 
487
      m2=istk(il2+1)
 
488
      top=top-1
 
489
c     get arg1
 
490
      il1=iadr(lstk(top))
 
491
      ilrs=il1
 
492
      if(istk(il1).lt.0) il1=iadr(istk(il1+1))
 
493
      m1=istk(il1+1)
 
494
 
 
495
      if(mn3.eq.0) then 
 
496
c     .  arg3=[]
 
497
         ilrs=iadr(lstk(top))
 
498
         istk(ilrs)=1
 
499
         istk(ilrs+1)=0
 
500
         istk(ilrs+2)=0
 
501
         istk(ilrs+3)=0
 
502
         lstk(top+1)=sadr(ilrs+4)+1
 
503
         goto 999
 
504
      elseif(m3.lt.0) then
 
505
c     .arg3=eye
 
506
         call error(14)
 
507
         return
 
508
      endif
 
509
c     check and convert indices variables
 
510
      call indxg(il1,m3,ili,mi,mxi,lw,1)
 
511
      if(err.gt.0) return
 
512
      if(mi.lt.0) then
 
513
         mr=mxi
 
514
      else
 
515
         mr=mi
 
516
      endif
 
517
      call indxg(il2,n3,ilj,nj,mxj,lw,1)
 
518
      if(err.gt.0) return
 
519
      if(nj.lt.0) then
 
520
         nr=mxj
 
521
      else
 
522
         nr=nj
 
523
      endif
 
524
c
 
525
 76   continue
 
526
      mn=mr*nr
 
527
      if(mn.eq.0) then 
 
528
c     .  arg1=[] or arg2=[] 
 
529
         ilrs=iadr(lstk(top))
 
530
         istk(ilrs)=1
 
531
         istk(ilrs+1)=0
 
532
         istk(ilrs+2)=0
 
533
         istk(ilrs+3)=0
 
534
         lstk(top+1)=sadr(ilrs+4)+1
 
535
         goto 999
 
536
      endif
 
537
      lptr=iadr(lw)
 
538
      irc=lptr+m3+1
 
539
      lw=sadr(irc+mr)
 
540
      nelr=iadr(lstk(bot))-iadr(lw)
 
541
      if(nelr.le.0) then
 
542
         err=lw-lstk(bot)
 
543
         call error(17)
 
544
         return
 
545
      endif
 
546
      lw=sadr(irc+mr+nelr)
 
547
      nel=nelr
 
548
      call lspe2(m3,n3,nel3,istk(irc3),istk(ili),mi,
 
549
     $        istk(ilj),nj,mr,nr,nelr,istk(irc),istk(lptr),ierr)
 
550
      ilrs=iadr(lstk(top))
 
551
      istk(ilrs)=6
 
552
      istk(ilrs+1)=mr
 
553
      istk(ilrs+2)=nr
 
554
      istk(ilrs+3)=0
 
555
      istk(ilrs+4)=nelr
 
556
      call icopy(mr+nelr,istk(irc),1,istk(ilrs+5),1)
 
557
      l1=sadr(ilrs+5+mr+nelr)
 
558
      lstk(top+1)=l1
 
559
      go to 999
 
560
c      
 
561
c     insert
 
562
 80   continue
 
563
      if(rhs.eq.4) goto 90
 
564
c     arg3(arg1)=arg2
 
565
c     get arg3      
 
566
      il3=iadr(lstk(top))
 
567
      if(istk(il3).lt.0) il3=iadr(istk(il3+1))
 
568
      m3=istk(il3+1)
 
569
      n3=istk(il3+2)
 
570
      it3=istk(il3+3)
 
571
      if(istk(il3).eq.6) then
 
572
         nel3=istk(il3+4)
 
573
         irc3=il3+5
 
574
         l3=sadr(irc3+m3+nel3)
 
575
      else
 
576
         top=top0
 
577
         fin=-fin
 
578
         return
 
579
      endif
 
580
      mn3=m3*n3
 
581
 
 
582
c     get arg2
 
583
      top=top-1
 
584
      il2=iadr(lstk(top))
 
585
      if(istk(il2).lt.0) il2=iadr(istk(il2+1))
 
586
      m2=istk(il2+1)
 
587
      n2=istk(il2+2)
 
588
      it2=istk(il2+3)
 
589
      if(istk(il2).eq.6) then
 
590
         nel2=istk(il2+4)
 
591
         irc2=il2+5
 
592
         l2=sadr(irc2+m2+nel2)
 
593
      elseif(istk(il2).eq.4) then
 
594
         l2=il2+3
 
595
         nel2=m2*n2
 
596
      elseif(istk(il2).eq.1.and.m2*n2.eq.0) then
 
597
         l2=sadr(il2+4)
 
598
         nel2=m2*n2
 
599
         istk(il2)=4
 
600
         l2=il2+3
 
601
      else
 
602
         top=top0
 
603
         fin=-fin
 
604
         return
 
605
      endif
 
606
      mn2=m2*n2
 
607
 
 
608
c     get arg1
 
609
      top=top-1
 
610
      il1=iadr(lstk(top))
 
611
      ilrs=il1
 
612
      if(istk(il1).lt.0) il1=iadr(istk(il1+1))
 
613
      m1=istk(il1+1)
 
614
      n1=istk(il1+2)
 
615
 
 
616
 
 
617
      if (m2.eq.0) then
 
618
c     .  arg3(arg1)=[] -->[]
 
619
         if(m1.eq.-1) then
 
620
c     .    arg3(:)=[] 
 
621
            ilrs=iadr(lstk(top))
 
622
            istk(ilrs)=1
 
623
            istk(ilrs+1)=0
 
624
            istk(ilrs+2)=0
 
625
            istk(ilrs+3)=0
 
626
            lstk(top+1)=sadr(ilrs+4)+1
 
627
            goto 999
 
628
         elseif(m1.eq.0) then
 
629
c     .     arg3([])=[]  --> arg3
 
630
            call icopy(5+m3+nel3,istk(il3),1,istk(ilrs),1)
 
631
            l=sadr(ilrs+5+m3+nel3)
 
632
            lstk(top+1)=l
 
633
            goto 999
 
634
         else
 
635
c     .     arg3(arg1)=[]
 
636
            if(istk(il1).eq.4.and.m3.eq.m1.and.n3.eq.n1) then
 
637
               if(.not.isany(il1)) then
 
638
c     .           arg3([])=[]  --> arg3
 
639
                  call icopy(5+m3+nel3,istk(il3),1,istk(ilrs),1)
 
640
                  l=sadr(ilrs+5+m3+nel3)
 
641
                  lstk(top+1)=l
 
642
                  goto 999
 
643
               endif
 
644
            endif
 
645
c     .     arg3(arg1)=[] -->arg3(compl(arg1),:)
 
646
            if(m3.gt.1.and.n3.gt.1) then
 
647
c     .        call macro coded op to reshape and insert
 
648
               top=top0
 
649
               fin=-fin
 
650
               return
 
651
            else
 
652
               call indxgc(il1,mn3,ilr,mi,mx,lw)
 
653
               if(err.gt.0) return
 
654
               l2=l3
 
655
               n2=n3
 
656
               m2=m3
 
657
               mn2=m2*n2
 
658
               it2=it3
 
659
               nel2=nel3
 
660
               irc2=irc3
 
661
c     .     call extraction
 
662
               goto 72
 
663
            endif
 
664
         endif
 
665
      elseif(m2.lt.0.or.m3.lt.0) then
 
666
c     .  arg3=eye,arg2=eye
 
667
         call error(14)
 
668
         return
 
669
      elseif(m1.lt.0) then
 
670
c     .  arg3(:)=arg2 reshape arg2 according to arg3
 
671
         if(mn2.ne.mn3) then
 
672
            call error(15)
 
673
            return
 
674
         endif
 
675
         if(m2.ne.m3) then
 
676
            top=top0
 
677
            fin=-fin
 
678
            return
 
679
         endif
 
680
         ilrs=iadr(lstk(top))
 
681
         istk(ilrs)=6
 
682
         istk(ilrs+1)=m3
 
683
         istk(ilrs+2)=n3
 
684
         call icopy(2+m2+nel2,istk(il2+3),1,istk(ilrs+3),1)
 
685
         l1=sadr(ilrs+5+m2+nel2)
 
686
         lstk(top+1)=l1
 
687
         return
 
688
      elseif(m3.gt.1.and.n3.gt.1) then
 
689
c     .  arg3(arg1)=arg2 with arg3 not a vector
 
690
         top=top0
 
691
         fin=-fin
 
692
         return
 
693
      endif
 
694
      call indxg(il1,mn3,ili,mi,mxi,lw,1)
 
695
      if(err.gt.0) return
 
696
      if(mi.eq.0) then
 
697
c     .  arg3([])=arg2
 
698
         call error(15)
 
699
         return
 
700
      endif
 
701
      if(mi.ne.mn2) then
 
702
         call error(15)
 
703
         return
 
704
      endif
 
705
c     
 
706
      if (n3.gt.1.and.m3.gt.1) then
 
707
c     .  arg3 is not a vector
 
708
         if(n2.gt.1.and.m2.gt.1) then
 
709
            call error(15)
 
710
            return
 
711
         endif
 
712
         if(mxi.gt.m3*n3) then
 
713
            call error(21)
 
714
            return
 
715
         endif
 
716
         mr=m3
 
717
         nr=n3
 
718
      elseif (n3.le.1.and.n2.le.1) then
 
719
c     .  arg3 and arg2 are  column vectors
 
720
         mr=max(m3,mxi)
 
721
         nr=max(n3,1)
 
722
      elseif (m3.le.1.and.m2.le.1) then
 
723
c     .  row vectors
 
724
         nr=max(n3,mxi)
 
725
         mr=max(m3,1)
 
726
      else
 
727
c     .  arg3 and arg2 dimensions dont agree
 
728
         call error(15)
 
729
         return
 
730
      endif
 
731
 
 
732
c     set output sizes
 
733
      if (m3 .gt. 1.or.m1.lt.0) then
 
734
c     .  column vector
 
735
         m=mi
 
736
         n=-1
 
737
         mr = mi
 
738
         nr = 1
 
739
      else
 
740
c     .  row vector
 
741
         m=-1
 
742
         n=mi
 
743
         nr = mi
 
744
         mr = 1
 
745
      endif
 
746
c     
 
747
      lptr=iadr(lw)
 
748
      irc=lptr+mr+1
 
749
      nelr=iadr(lstk(bot))-irc-mr
 
750
      if(nelr.le.0) then
 
751
         err=lw-lstk(bot)
 
752
         call error(17)
 
753
         return
 
754
      endif
 
755
      lw=sadr(irc+mr+nelr)
 
756
      nel=nelr
 
757
      if(istk(il3).eq.6) then
 
758
         if(istk(il2).eq.6) then
 
759
            call lspisp(m3,n3,nel3,istk(irc3),istk(ili),m,istk(ili),n,m2
 
760
     $           ,n2,nel2,istk(irc2),mr,nr,nelr,istk(irc),istk(lptr)
 
761
     $           ,ierr)
 
762
         elseif(istk(il2).eq.4) then
 
763
            call lspis(m3,n3,nel3,istk(irc3),istk(ili),m,istk(ili),n,m2
 
764
     $           ,n2,istk(l2),mr,nr,nelr,istk(irc),ierr) 
 
765
         endif
 
766
      endif
 
767
      if(ierr.ne.0) then
 
768
         buf='not enough memory'
 
769
         call error(9999)
 
770
         return
 
771
      endif
 
772
      ilrs=iadr(lstk(top))
 
773
      istk(ilrs)=6
 
774
      istk(ilrs+1)=mr
 
775
      istk(ilrs+2)=nr
 
776
      istk(ilrs+3)=0
 
777
      istk(ilrs+4)=nelr
 
778
      call icopy(mr+nelr,istk(irc),1,istk(ilrs+5),1)
 
779
      l1=sadr(ilrs+5+mr+nelr)
 
780
      lstk(top+1)=l1
 
781
      go to 999
 
782
c     
 
783
 90   continue
 
784
c     arg4(arg1,arg2)=arg3
 
785
c     get arg4      
 
786
      il4=iadr(lstk(top))
 
787
      if(istk(il4).lt.0) il4=iadr(istk(il4+1))
 
788
      m4=istk(il4+1)
 
789
      n4=istk(il4+2)
 
790
      it4=istk(il4+3)
 
791
      if(istk(il4).eq.6) then
 
792
         nel4=istk(il4+4)
 
793
         irc4=il4+5
 
794
         l4=sadr(irc4+m4+nel4)
 
795
      else
 
796
         top=top0
 
797
         fin=-fin
 
798
         return
 
799
      endif
 
800
      mn4=m4*n4
 
801
c     get arg3
 
802
      top=top-1
 
803
      il3=iadr(lstk(top))
 
804
      if(istk(il3).lt.0) il3=iadr(istk(il3+1))
 
805
      m3=istk(il3+1)
 
806
      n3=istk(il3+2)
 
807
      it3=istk(il3+3)
 
808
      if(istk(il3).eq.6) then
 
809
         nel3=istk(il3+4)
 
810
         irc3=il3+5
 
811
         l3=sadr(irc3+m3+nel3)
 
812
      elseif(istk(il3).eq.4) then
 
813
         nel3=m3*n3
 
814
         l3=il3+3
 
815
      elseif(istk(il3).eq.1.and.m3*n3.eq.0) then
 
816
         l3=sadr(il3+4)
 
817
         nel3=m3*n3
 
818
         istk(il3)=4
 
819
         l3=il3+3
 
820
      else
 
821
         top=top0
 
822
         fin=-fin
 
823
         return
 
824
      endif
 
825
      mn3=m3*n3
 
826
c     get arg2
 
827
      top=top-1
 
828
      il2=iadr(lstk(top))
 
829
      if(istk(il2).lt.0) il2=iadr(istk(il2+1))
 
830
      m2=istk(il2+1)
 
831
c     get arg1
 
832
      top=top-1
 
833
      il1=iadr(lstk(top))
 
834
      ilrs=il1
 
835
      if(istk(il1).lt.0) il1=iadr(istk(il1+1))
 
836
      m1=istk(il1+1)
 
837
 
 
838
      if (m3.eq.0) then
 
839
c     .  arg4(arg1,arg2)=[]
 
840
         if(m1.eq.-1.and.m2.eq.-1) then
 
841
c     .    arg4(:,:)=[] -->[]
 
842
            ilrs=iadr(lstk(top))
 
843
            istk(ilrs)=1
 
844
            istk(ilrs+1)=0
 
845
            istk(ilrs+2)=0
 
846
            istk(ilrs+3)=0
 
847
            lstk(top+1)=sadr(ilrs+4)+1
 
848
            goto 999
 
849
         elseif(m1.eq.0.or.m2.eq.0) then
 
850
c     .     arg4([],arg2)=[],  arg4(arg1,[])=[] --> arg4
 
851
            call icopy(5+m4+nel4,istk(il4),1,istk(ilrs),1)
 
852
            l=sadr(ilrs+5+m4+nel4)
 
853
            lstk(top+1)=l
 
854
            goto 999
 
855
         elseif(m2.eq.-1) then
 
856
c     .     arg3(arg1,:)=[] --> arg3(compl(arg1),:)
 
857
            call indxgc(il1,m4,ili,mi,mxi,lw)
 
858
            if(err.gt.0) return
 
859
            mr=mi
 
860
            call indxg(il2,n4,ilj,nj,mxj,lw,1)
 
861
            if(err.gt.0) return
 
862
            if(nj.lt.0) then
 
863
               nr=mxj
 
864
            else
 
865
               nr=nj
 
866
            endif
 
867
            l3=l4
 
868
            n3=n4
 
869
            m3=m4
 
870
            mn3=m3*n3
 
871
            it3=it4
 
872
            irc3=irc4
 
873
            nel3=nel4
 
874
c     .     call extraction
 
875
            goto 76
 
876
         elseif(m1.eq.-1) then
 
877
c     .     arg3(:,arg2)=[] --> arg3(:,compl(arg2))
 
878
            call indxgc(il2,n4,ilj,nj,mxj,lw)
 
879
            if(err.gt.0) return
 
880
            nr=nj
 
881
            call indxg(il1,m4,ili,mi,mxi,lw,1)
 
882
            if(err.gt.0) return
 
883
            if(mi.lt.0) then
 
884
               mr=mxi
 
885
            else
 
886
               mr=mi
 
887
            endif
 
888
            l3=l4
 
889
            n3=n4
 
890
            m3=m4
 
891
            mn3=m3*n3
 
892
            it3=it4
 
893
            irc3=irc4
 
894
            nel3=nel4
 
895
c     .     call extraction
 
896
            goto 76
 
897
         else
 
898
c     .     arg4(arg1,arg2)=[] 
 
899
            lw1=lw
 
900
            call indxgc(il2,n4,ilj,nj,mxj,lw)
 
901
            if(err.gt.0) return
 
902
            nr=nj
 
903
            if(nj.eq.0) then
 
904
c     .        arg4(arg1,1:n4)=[] 
 
905
               call indxgc(il1,m4,ili,mi,mxi,lw)
 
906
               lw2=lw
 
907
               if(err.gt.0) return
 
908
               mr=mi
 
909
c     .        arg2=1:n4
 
910
               if(mi.eq.0) then
 
911
c     .           arg4(1:m4,1:n4)=[] 
 
912
                  ilrs=iadr(lstk(top))
 
913
                  istk(ilrs)=1
 
914
                  istk(ilrs+1)=0
 
915
                  istk(ilrs+2)=0
 
916
                  istk(ilrs+3)=0
 
917
                  lstk(top+1)=sadr(ilrs+4)+1
 
918
                  goto 999
 
919
               else
 
920
c     .           arg4(arg1,1:n4)=[] 
 
921
c     .           replace arg2 by ":"
 
922
                  il2=iadr(lw2)
 
923
                  istk(il2)=1
 
924
                  istk(il2+1)=-1
 
925
                  istk(il2+2)=-1
 
926
                  istk(il2+3)=0
 
927
c     .
 
928
                  lw=lw2+2
 
929
                  call indxg(il2,n4,ilj,nj,mxj,lw,1)
 
930
                  if(err.gt.0) return
 
931
                  if(nj.lt.0) then
 
932
                     nr=mxj
 
933
                  else
 
934
                     nr=nj
 
935
                  endif
 
936
                  l3=l4
 
937
                  n3=n4
 
938
                  m3=m4
 
939
                  it3=it4
 
940
                  mn3=m3*n3
 
941
                  irc3=irc4
 
942
                  nel3=nel4
 
943
c     .           call extraction
 
944
                  goto 76
 
945
               endif
 
946
            else
 
947
               lw=lw1
 
948
               call indxgc(il1,m4,ili,mi,mxi,lw)
 
949
               if(err.gt.0) return
 
950
               if(mi.eq.0) then
 
951
c     .           arg4(1:m4,arg2)=[] 
 
952
                  call indxg(il1,m4,ili,mi,mxi,lw,1)
 
953
                  if(err.gt.0) return
 
954
                  if(mi.lt.0) then
 
955
                     mr=mxi
 
956
                  else
 
957
                     mr=mi
 
958
                  endif
 
959
                  l3=l4
 
960
                  n3=n4
 
961
                  m3=m4
 
962
                  it3=it4
 
963
                  mn3=m3*n3
 
964
                  irc3=irc4
 
965
                  nel3=nel4
 
966
c     .           call extraction
 
967
                  goto 76
 
968
               else
 
969
                  call error(15)
 
970
                  return
 
971
               endif
 
972
            endif
 
973
         endif
 
974
      elseif(m3.lt.0.or.m4.lt.0) then
 
975
c     .  arg3=eye , arg4=eye
 
976
         call error(14)
 
977
         return
 
978
      elseif(m1.eq.-1.and.m2.eq.-1) then
 
979
c     .  arg4(:,:)=arg3
 
980
         if(mn3.ne.mn4) then
 
981
            call error(15)
 
982
            return
 
983
         endif
 
984
         if(m3.ne.m4) then
 
985
            top=top0
 
986
            fin=-fin
 
987
            return
 
988
         endif
 
989
c     .  reshape arg3 according to arg4
 
990
         ilrs=iadr(lstk(top))
 
991
         istk(ilrs)=6
 
992
         istk(ilrs+1)=m4
 
993
         istk(ilrs+2)=n4
 
994
         call icopy(2+m3+nel3,istk(il3+3),1,istk(ilrs+3),1)
 
995
         l1=sadr(ilrs+5+m3+nel3)
 
996
         lstk(top+1)=l1
 
997
         return
 
998
      endif
 
999
 
 
1000
      call indxg(il1,m4,ili,mi,mxi,lw,1)
 
1001
      if(err.gt.0) return
 
1002
      if(mi.lt.0) then
 
1003
         mr1=mxi
 
1004
      else
 
1005
         mr1=mi
 
1006
      endif
 
1007
      call indxg(il2,n4,ilj,mj,mxj,lw,1)
 
1008
      if(err.gt.0) return
 
1009
      if(mj.lt.0) then
 
1010
         nr1=mxj
 
1011
      else
 
1012
         nr1=mj
 
1013
      endif
 
1014
      if(mr1.eq.0.or.nr1.eq.0) then
 
1015
         call error(15)
 
1016
         return
 
1017
      endif
 
1018
      if(mr1.ne.m3.or.nr1.ne.n3) then
 
1019
c     .  sizes of arg1 or arg2 dont agree with arg3 sizes
 
1020
         call error(15)
 
1021
         return
 
1022
      endif
 
1023
      mr=max(m4,mxi)
 
1024
      nr=max(n4,mxj)
 
1025
c
 
1026
      lptr=iadr(lw)
 
1027
      irc=lptr+mr+1
 
1028
      nelr=iadr(lstk(bot))-irc-mr
 
1029
      if(nelr.le.0) then
 
1030
         err=lw-lstk(bot)
 
1031
         call error(17)
 
1032
         return
 
1033
      endif
 
1034
      lw=sadr(irc+mr+nelr)
 
1035
      nel=nelr
 
1036
      if(istk(il4).eq.6) then
 
1037
         if(istk(il3).eq.6) then
 
1038
            call lspisp(m4,n4,nel4,istk(irc4),istk(ili),mi,istk(ilj),mj,
 
1039
     $           m3,n3,nel3,istk(irc3),mr,nr,nelr,istk(irc),istk(lptr)
 
1040
     $           ,ierr)
 
1041
         elseif(istk(il3).eq.4) then
 
1042
            l3=il3+3
 
1043
            call lspis(m4,n4,nel4,istk(irc4),istk(ili),mi,istk(ilj),mj
 
1044
     $           ,m3,n3,istk(l3),mr,nr,nelr,istk(irc),ierr) 
 
1045
         endif
 
1046
      endif
 
1047
      if(ierr.ne.0) then
 
1048
         buf='not enough memory'
 
1049
         call error(9999)
 
1050
         return
 
1051
      endif
 
1052
      ilrs=iadr(lstk(top))
 
1053
      istk(ilrs)=6
 
1054
      istk(ilrs+1)=mr
 
1055
      istk(ilrs+2)=nr
 
1056
      istk(ilrs+3)=0
 
1057
      istk(ilrs+4)=nelr
 
1058
      call icopy(mr+nelr,istk(irc),1,istk(ilrs+5),1)
 
1059
      l1=sadr(ilrs+5+mr+nelr)
 
1060
      lstk(top+1)=l1
 
1061
      go to 999
 
1062
c
 
1063
 130  continue
 
1064
c     comparaisons
 
1065
      if(mn2.eq.0.and.mn1.eq.0) then
 
1066
         if(op.eq.equal.or.op.eq.less+great) then
 
1067
            il1=iadr(lstk(top))
 
1068
            istk(il1)=4
 
1069
            istk(il1+1)=1
 
1070
            istk(il1+2)=1
 
1071
            istk(il1+3)=1
 
1072
            if(op.eq.less+great) istk(il1+3)=0
 
1073
            lstk(top+1)=sadr(il1+4)
 
1074
            goto 999
 
1075
         else
 
1076
            call error(60)
 
1077
            return
 
1078
         endif
 
1079
      endif
 
1080
      if(mn1.ne.1.and.mn2.ne.1) then
 
1081
         if(n1.ne.n2.or.m1.ne.m2) then
 
1082
            if(op.eq.equal.or.op.eq.less+great) then
 
1083
               il1=iadr(lstk(top))
 
1084
               istk(il1)=4
 
1085
               istk(il1+1)=1
 
1086
               istk(il1+2)=1
 
1087
               istk(il1+3)=0
 
1088
               if(op.eq.less+great) istk(il1+3)=1
 
1089
               lstk(top+1)=sadr(il1+4)
 
1090
               return
 
1091
            else
 
1092
               call error(60)
 
1093
               return
 
1094
            endif
 
1095
         endif
 
1096
      endif
 
1097
      if(istk(il1).ne.4.and.istk(il1).ne.6.or.istk(il2).ne.4.and
 
1098
     $     .istk(il2).ne.6) then
 
1099
         top=top0
 
1100
         fin=-fin
 
1101
         return
 
1102
      endif
 
1103
c
 
1104
      mr=m1
 
1105
      nr=n1
 
1106
      if(m1*n1.eq.1) then
 
1107
         mr=m2
 
1108
         nr=n2
 
1109
      endif
 
1110
      irc=iadr(lw)
 
1111
      nelmx=(iadr(lstk(bot))-irc-mr-10)
 
1112
      lw=sadr(irc+mr+nelmx)
 
1113
      err=lw-lstk(bot)
 
1114
      if(err.gt.0) then
 
1115
         call error(17)
 
1116
         return
 
1117
      endif  
 
1118
      nel=nelmx
 
1119
      if(istk(il1).eq.4) then
 
1120
         l1=il1+3
 
1121
         call lsosp(op,m1,n1,istk(l1),m2,n2,nel2,istk(irc2),nel,istk(irc
 
1122
     $        ),ierr)
 
1123
      elseif(istk(il2).eq.4) then
 
1124
         l2=il2+3
 
1125
         call lspos(op,m1,n1,nel1,istk(irc1),
 
1126
     $           m2,n2,istk(l2),nel,istk(irc),ierr)
 
1127
      else
 
1128
         call lsposp(op,m1,n1,nel1,istk(irc1),
 
1129
     $           m2,n2,nel2,istk(irc2),
 
1130
     $           nel,istk(irc),ierr)
 
1131
      endif
 
1132
      if(ierr.ne.0) then
 
1133
         buf='not enough memory'
 
1134
         call error(9999)
 
1135
         return
 
1136
      endif
 
1137
      il1=iadr(lstk(top))
 
1138
      istk(il1)=6
 
1139
      istk(il1+1)=mr
 
1140
      istk(il1+2)=nr
 
1141
      istk(il1+3)=0
 
1142
      istk(il1+4)=nel
 
1143
      irc1=il1+5
 
1144
      call icopy(mr+nel,istk(irc),1,istk(irc1),1)
 
1145
      l1=sadr(irc1+mr+nel)
 
1146
      lstk(top+1)=l1
 
1147
      go to 999
 
1148
c     
 
1149
 999  return
 
1150
      end
 
1151
 
 
1152
 
 
1153