From bug-octave-request  Mon Dec 21 07:46:40 1992
Subject: core dump
From: jbraw (James B. Rawlings)
To: bug-octave
Date: Mon, 21 Dec 92 07:46:38 -0600


John, I finally got one.  Here is the history minus the last line that
caused the core dump.  The last line was (to the best of my knowledge)

ans=x*c; save x.tmp ans; ans=xhat*c; save xhat.tmp ans;

I can give you a briefer list of the history that should give you the
same error.  I appended dlkqe.m and simdlkqe.m, they are pretty short.

--octav_hist--
his
history
quit
help history
history
a=zeros(10,10); a(6:10,6:10)=eye(5); c=[eye(5) eye(5);g=eye(10);
eye(5)
a=zeros(10,10); a(6:10,6:10)=eye(5); c=[eye(5) eye(5)];g=eye(10);
q=eye(10);r=1;niter=10;p0=eye(10);icp=[1 2 3 4 5 4 3 2]';
history -w
[l,p]=dlkqe(a,c,g,q,r,icp,p0,niter)
l
icp=[1 2 3 4 5 5 4 3 2 1]'
[l,p]=dlkqe(a,c,g,q,r,icp,p0,niter)
l
save l.tmp l; save p.tmp p
[l,p]=dlkqe(a,c,g,q,r,icp,p0,niter);
l
p
icp=[1 2 3 4 5 4 3 2]'
[l,p]=dlkqe(a,c,g,q,r,icp,p0,niter);
l
p
g
g(1:5,1:5)=zeros(5)
[l,p]=dlkqe(a,c,g,q,r,icp,p0,niter);
l
icp=[1 1 1 1 1]';
[l,p]=dlkqe(a,c,g,q,r,icp,p0,niter);
l
g=eye(10);icp=[1 2 3 4 5 4 3 2]';
history -w
history
x0=zeros(1,10);x=zeros(1,10);ntimes=5;
[x,xhat,y]=simdlkqe(a,c,g,l,x0,xhat0,ntimes);
x0=zeros(1,10);xhat0=zeros(1,10);ntimes=5;
[x,xhat,y]=simdlkqe(a,c,g,l,x0,xhat0,ntimes);
l
icp
[l,p]=dlkqe(a,c,g,q,r,icp,p0,niter);
l
[l,p]=dlkqe(a,c,g,q,r,icp,p0,niter);
[x,xhat,y]=simdlkqe(a,c,g,l,x0,xhat0,ntimes);
c
[x,xhat,y]=simdlkqe(a,c,g,l,icp,x0,xhat0,ntimes);
x
rand('normal');
[x,xhat,y]=simdlkqe(a,c,g,l,icp,x0,xhat0,ntimes);
x
y
xhat
save x.tmp x; save xhat.tmp  xhat;
history -w

--dlkqe.m--
function [l,p] = dlkqe(a,c,g,q,r,icp,p0,niter)
%
%  Warning: this code is not debugged, has absolutely no error checking
%  or comments.  Do not use this is if your name is not:
%
%  James B. Rawlings
%  Department of Chemical Engineering
%  University of Texas at Austin
%  Austin, Texas 78712
%
%  December 18, 1992
%
%  if you do, you will regret it.
%
[n n]=size(a);
%
%  nsensors=1 is the current set up.  Can add sensors later if necessary.
%
[nperiod nsensors]=size(icp);
pk=p0;
for i=1:niter
for k=1:nperiod
%
%  scanning sensor, update which locations on the sheet are being 
%  sensed at the kth time through the cycle.
%
  ck=zeros(1,n/2);
  ck(1,icp(k)) = 1;
  ck=ck*c;
  pk1 = a*pk*a' - a*pk*ck'*inv(ck*pk*ck'+r)*ck*pk*a' + g*q*g';
  p( (k-1)*n+1:k*n, 1:n ) = pk1;
  l( k, 1:n ) = (a*pk1*ck'*inv(ck*pk1*ck'+r))';
  pk=pk1;
endfor
endfor
endfunction

--simdlkqe.m---
function [x,xhat,y] = simdlkqe(a,c,g,l,icp,x0,xhat0,ntimes)
%
%  Warning: this code is not debugged, has absolutely no error checking
%  or comments.  Do not use this is if your name is not:
%
%  James B. Rawlings
%  Department of Chemical Engineering
%  University of Texas at Austin
%  Austin, Texas 78712
%
%  December 18, 1992
%
%  if you do, you will regret it.
%
[n n]=size(a);
[nperiod nsensors]=size(icp);
index=1;
x(index,1:n)=x0;
xhat(index,1:n)=xhat0;
for i=1:ntimes
for k=1:nperiod
%
%  scanning sensor, update which locations on the sheet are being 
%  sensed at the kth time through the cycle.
%
  ck=zeros(1,n/2);
  ck(1,icp(k)) = 1;
  ck=ck*c;
  lk = l( k, 1:n )';
  y(index) = ck*x(index,1:n)' + rand;
  x(index+1,1:n) = x(index,1:n)*a' + rand(1,n)*g';
  ek = y(index) - ck*xhat(index,1:n)';
  xhat(index+1,1:n) = xhat(index,1:n)*a' + ek'*lk';
  index = index + 1;
endfor
endfor
endfunction