From bug-octave-request at che dot utexas dot edu  Fri Jun  3 15:56:52 1994
Subject: Weird problem
From: charles dot henkel at srs dot gov (Chuck Henkel)
To: bug-octave at che dot utexas dot edu
Date: Fri, 03 Jun 1994 16:55:45 -0400

Context:  Octave 1.0
	  RS/6000 AIX 3.2.5
	  SGI Irix 5.2 (Running 4.0.5 binary)

Problem: 

  Here is a small m-file. If "sgn" is 1., then it works OK. If "sgn"
  is -1., then, after a minute or so of crunching, I get the message
  "Virtual memory exceeded in `new'" and octave crashes.

  Note that sgn is simply used to decide which root of a quadratic I
  choose.

sgn=-1.;
nit=100;
n=10;
nn=n*n;
siga=0.5;
sigf=0.8;
A=mkmat(n,siga);
F=sigf*diag(ones(nn,1));
R=diag(diag(A)) + diag(diag(A,-1),-1) + diag(diag(A,1),1);

x=ones(nn,1);
xAx=x'*A*x;
xFx=x'*F*x;
sig=xAx/xFx;

for i = 2:nit
  r=(sig*F-A)*x;
  g=-(2.*r)/xFx;
  z=R\r;
  gtz = g'*z;
  if i > 2 
    beta = gtz/gtzo;
    p = -z + beta*p;
  else
    gtzo = gtz;
    p = -z;
  end
  pAp = p'*A*p;
  xAp = x'*A*p;
  pFp = p'*F*p;
  xFp = x'*F*p;
  u = pAp*xFp - xAp*pFp;
  v = pAp*xFx - xAx*pFp;
  w = xAp*xFx - xAx*xFp;
  alpha = (-v + sgn*sqrt(v*v-4.*u*w))/(2.*u);
  x = x + alpha*p;
  xAx=x'*A*x;
  xFx=x'*F*x;
  sig=xAx/xFx
end

%
% This is the auxiliary function mkmat
%

function a = mkmat(n,sigma)
% Generate a five banded matrix of dimension n*n
% Zero flux bc on W, S
% Zero current on E, N
nn = n*n;
[ic,ir] = size(sigma);
irc = ic*ir;

if (irc == 1),
 for k=1:nn
       a(k,k) = 4. + sigma;
 end
elseif (irc == nn),
 for k=1:nn
       a(k,k) = 4. + sigma(k);
 end
end
for k=1:nn-1
      a(k,k+1) = -1.;
      a(k+1,k) = -1.;
end
for k=n+1:n:nn
      a(k,k-1) = 0.;
end
for k=n:n:nn-n
      a(k,k+1) = 0.;
end
for k=1:nn-n
      a(k,k+n) = -1.;
      a(k+n,k) = -1.;
end
for k=1:n:nn
      a(k,k) = a(k,k) + 1;
end
for k=n:n:nn
      a(k,k) = a(k,k) - 1;
end
for k=1:n
      a(k,k) = a(k,k) + 1;
end
for k=nn-n+1:nn
      a(k,k) = a(k,k) - 1;
end