From bug-octave-request at che dot utexas dot edu  Sat Sep  3 17:42:49 1994
Subject: Another crash
From: Edward Scott Meadows <esm at che dot utexas dot edu>
To: bug-octave
Date: Sat, 3 Sep 1994 17:42:40 -0500

I can't tell what happened here.  It seems to be a crash in NPSOL.
The function minconst is in /home/esm/octave/k1k2/ and given below.
Sorry, but I lost A when octave crashed and couldn't recover it from
the screen buffer.  It was created using A=rand(5)*4-2.

octave:117> d=rand(4,1)
d =

  0.59189
  0.20789
  0.12786
  0.83189

octave:118> D=rand(4,5)*2-1
D =

  -0.491602   0.451000  -0.203817   0.669610   0.714030
   0.619066  -0.408047  -0.979774  -0.902463   0.054560
  -0.470024   0.847183  -0.694179   0.825065   0.061696
   0.768748  -0.264256  -0.987411  -0.787296   0.132401

octave:119> minconst (A,D,d)
forrtl: severe: recursive I/O operation, unit 6, file (null)
june:95> 

---------------------Here is minconst------------------------------------

function [C,N]=minconst(A,D,d);

# usage: [C,N]=minconst(A,D,d)
#
#    For the system
#
#    	x_{k+1} = A x_k
#    
#    subject to constraints
#
#    	D x_k <= d    k = 0,1,2 ... 
#
#    this function provides the minimal N for which the
#    constraints are guaranteed to be satisfied for all k >= N
#    and provides the corresponding system of constraints
#
#    	C = [ D ; D A ; D A^2 ; ... D A^N ]
#
#    such that C x_0 <= d guarantees D x_k <= d for all k. 
#    This is useful for linear systems to determine whether an 
#    initial condition will remain in a linearly constrained
#    region as the time index k increases.
#
#    Algorithm is from Gilbert and Tan, IEEE-TAC(36):9,1991.
#    Users must beware that finite N cannot be expected for 
#    arbitrary (A,D,d).

[NA,MA]=size(A);
[ND,MD]=size(D);

# Error checking on matrix dimensions

if (NA!=MA)
  printf("minconst: Matrix in first argument non-square \n");
  C=[];
  N=0;
  return;
endif;

if (MD!=NA||ND!=length(d))
  printf("minconst: Inconsistent dimensions on input \n");
  C=[];
  N=0;
  return;
endif;

# Algorithm requires solution of linear program.
# Sets up input for LP and begins to assemble output variables.

Alp=D;
ubl=d;
lbl=-Inf*ones(length(d),1);

ub=Inf*ones(MA,1);
lb=-Inf*ones(MA,1);

N=0;
Clp=D*A;
u=0*ones(NA,1);
H=0*ones(NA);
page_screen_output='false';
istop=1;

# Executes recursive series of linear programs to implement
# Gilbert and Tan algorithm.  (Have to use QPSOL since LPSOL
# not implemented in octave yet.)

while (istop==1)
  maxvec=0*ones(ND,1);
  for j=1:ND;
     [u,obj,info]=qpsol(u,H,-Clp(j,:),lb,ub,lbl,Alp,ubl);
     maxvec(j)=-obj-d(j);
  end;

  sum((maxvec>0));

# End condition is when max objective value below bound.  
# This is tested by next line.  I also insist that QPSOL
# exit properly with info = 0 or 1.

  if (sum((maxvec>0))==0&&info<2) 
     istop=0;

# This assembles the matrices for the next recursion.

   else
     Alp=[Alp;Clp];
     Clp=Clp*A;
     lbl=[lbl;-Inf*d];
     ubl=[ubl;d];
     N=N+1;
  endif;
endwhile;

page_screen_output='true';
C=Alp;
return;
end;