[m,n,a,ia,niv,iflag,s,p,ip,inhomo,v,w,wcnd]=mgsbv(m,n,a,ia,niv,iflag,s,p,ip,inhomo,v,w,wcnd);

function [m,n,a,ia,niv,iflag,s,p,ip,inhomo,v,w,wcnd]=mgsbv(m,n,a,ia,niv,iflag,s,p,ip,inhomo,v,w,wcnd);
persistent dot i ip1 ix iz j jk jp jq jy jz k kd kj kp l lix lr m2 nivn nmnr nn np1 nr nrm1 pjp psave ry sv t vl vnorm y ; 

global ml18jr_1; if isempty(ml18jr_1), ml18jr_1=0; end;
if isempty(dot), dot=0; end;
global ml5mco_3; if isempty(ml5mco_3), ml5mco_3=0; end;
global ml5mco_6; if isempty(ml5mco_6), ml5mco_6=0; end;
if isempty(pjp), pjp=0; end;
if isempty(psave), psave=0; end;
global ml18jr_2; if isempty(ml18jr_2), ml18jr_2=0; end;
if isempty(ry), ry=0; end;
global ml5mco_4; if isempty(ml5mco_4), ml5mco_4=0; end;
global ml5mco_2; if isempty(ml5mco_2), ml5mco_2=0; end;
if isempty(sv), sv=0; end;
if isempty(t), t=0; end;
global ml18jr_3; if isempty(ml18jr_3), ml18jr_3=0; end;
global ml5mco_5; if isempty(ml5mco_5), ml5mco_5=0; end;
global ml5mco_1; if isempty(ml5mco_1), ml5mco_1=0; end;
if isempty(vl), vl=0; end;
if isempty(vnorm), vnorm=0; end;
if isempty(y), y=0; end;
if isempty(i), i=0; end;
global ml18jr_18; if isempty(ml18jr_18), ml18jr_18=0; end;
global ml18jr_11; if isempty(ml18jr_11), ml18jr_11=0; end;
global ml18jr_12; if isempty(ml18jr_12), ml18jr_12=0; end;
if isempty(ip1), ip1=0; end;
if isempty(ix), ix=0; end;
if isempty(iz), iz=0; end;
if isempty(j), j=0; end;
if isempty(jk), jk=0; end;
if isempty(jp), jp=0; end;
if isempty(jq), jq=0; end;
if isempty(jy), jy=0; end;
if isempty(jz), jz=0; end;
if isempty(k), k=0; end;
if isempty(kd), kd=0; end;
if isempty(kj), kj=0; end;
if isempty(kp), kp=0; end;
if isempty(l), l=0; end;
if isempty(lix), lix=0; end;
global ml5mco_7; if isempty(ml5mco_7), ml5mco_7=0; end;
if isempty(lr), lr=0; end;
if isempty(m2), m2=0; end;
global ml18jr_7; if isempty(ml18jr_7), ml18jr_7=0; end;
global ml18jr_8; if isempty(ml18jr_8), ml18jr_8=0; end;
global ml18jr_10; if isempty(ml18jr_10), ml18jr_10=0; end;
global ml18jr_15; if isempty(ml18jr_15), ml18jr_15=0; end;
global ml18jr_17; if isempty(ml18jr_17), ml18jr_17=0; end;
global ml18jr_5; if isempty(ml18jr_5), ml18jr_5=0; end;
if isempty(nivn), nivn=0; end;
if isempty(nmnr), nmnr=0; end;
if isempty(nn), nn=0; end;
global ml18jr_6; if isempty(ml18jr_6), ml18jr_6=0; end;
if isempty(np1), np1=0; end;
global ml18jr_13; if isempty(ml18jr_13), ml18jr_13=0; end;
if isempty(nr), nr=0; end;
if isempty(nrm1), nrm1=0; end;
global ml18jr_9; if isempty(ml18jr_9), ml18jr_9=0; end;
global ml18jr_14; if isempty(ml18jr_14), ml18jr_14=0; end;
global ml18jr_16; if isempty(ml18jr_16), ml18jr_16=0; end;
global ml18jr_4; if isempty(ml18jr_4), ml18jr_4=0; end;
%***BEGIN PROLOGUE  MGSBV
%***SUBSIDIARY
%***PURPOSE  Subsidiary to BVSUP
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (MGSBV-S, DMGSBV-D)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
% **********************************************************************
% Orthogonalize a set of N real vectors and determine their rank
%
% **********************************************************************
% INPUT
% **********************************************************************
%   M = Dimension of vectors
%   N = No. of vectors
%   A = Array whose first N cols contain the vectors
%   IA = First dimension of array A (col length)
%   NIV = Number of independent vectors needed
%   INHOMO = 1 Corresponds to having a non-zero particular solution
%   V = Particular solution vector (not included in the pivoting)
%   INDPVT = 1 Means pivoting will not be used
%
% **********************************************************************
% OUTPUT
% **********************************************************************
%   NIV = No. of linear independent vectors in input set
%     A = Matrix whose first NIV cols. contain NIV orthogonal vectors
%         which span the vector space determined by the input vectors
%   IFLAG
%          = 0 success
%          = 1 incorrect input
%          = 2 rank of new vectors less than N
%   P = Decomposition matrix.  P is upper triangular and
%             (old vectors) = (new vectors) * P.
%         The old vectors will be reordered due to pivoting
%         The dimension of p must be .GE. N*(N+1)/2.
%             (  N*(2*N+1) when N .NE. NFCC )
%   IP = Pivoting vector. The dimension of IP must be .GE. N.
%             (  2*N when N .NE. NFCC )
%   S = Square of norms of incoming vectors
%   V = Vector which is orthogonal to the vectors of A
%   W = Orthogonalization information for the vector V
%   WCND = Worst case (smallest) norm decrement value of the
%          vectors being orthogonalized  (represents a test
%          for linear dependence of the vectors)
% **********************************************************************
%
%***SEE ALSO  BVSUP
%***ROUTINES CALLED  PRVEC, SDOT
%***COMMON BLOCKS    ML18JR, ML5MCO
%***REVISION HISTORY  (YYMMDD)
%   750601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   890921  Realigned order of variables in certain COMMON blocks.
%           (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%   910722  Updated AUTHOR section.  (ALS)
%***end PROLOGUE  MGSBV
%
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([ia])).*prod([ia])-numel(a))],ia,[]);
v_shape=size(v);v=reshape(v,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
p_shape=size(p);p=reshape(p,1,[]);
ip_shape=size(ip);ip=reshape(ip,1,[]);
s_shape=size(s);s=reshape(s,1,[]);
%
%
% common :: ;
%% common /ml18jr/ ae , re , tol , nxpts , nic , nopg , mxnon ,ndisk , ntape , neq , indpvt , integ , nps , ntp ,neqivp , numort , nfcc , icoco;
%% common /ml18jr/ ml18jr_1 , ml18jr_2 , ml18jr_3 , ml18jr_4 , ml18jr_5 , ml18jr_6 , ml18jr_7 ,ml18jr_8 , ml18jr_9 , ml18jr_10 , ml18jr_11 , ml18jr_12 , ml18jr_13 , ml18jr_14 ,ml18jr_15 , ml18jr_16 , ml18jr_17 , ml18jr_18;
%
% common :: ;
%% common /ml5mco/ uro , sru , eps , sqovfl , twou , fouru , lpar;
%% common /ml5mco/ ml5mco_1 , ml5mco_2 , ml5mco_3 , ml5mco_4 , ml5mco_5 , ml5mco_6 , ml5mco_7;
%
%***FIRST EXECUTABLE STATEMENT  MGSBV
if( m>0 && n>0 && ia>=m )
%
jp = 0;
iflag = 0;
np1 = fix(n + 1);
y = 0.0;
m2 = fix(fix(m./2));
%
%     CALCULATE SQUARE OF NORMS OF INCOMING VECTORS AND SEARCH FOR
%     VECTOR WITH LARGEST MAGNITUDE
%
j = 0;
for i = 1 : n;
[vl ,m,dumvar2,dumvar4,dumvar2]=sdot(m,a(sub2ind(size(a),1,i):end),1,a(sub2ind(size(a),1,i):end),1);      a(sub2ind(size(a),1,i):end)=dumvar2; 
s(i) = vl;
if( n~=ml18jr_17 )
j = fix(2.*i - 1);
p(j) = vl;
ip(j) = fix(j);
end;
j = fix(j + 1);
p(j) = vl;
ip(j) = fix(j);
if( vl>y )
y = vl;
ix = fix(i);
end;
end; i = fix(n+1);
if( ml18jr_11==1 )
ix = 1;
y = p(1);
end;
lix = fix(ix);
if( n~=ml18jr_17 )
lix = fix(2.*ix - 1);
end;
p(lix) = p(1);
s(np1) = 0.;
if( inhomo==1 )
v_orig=v;    [ s(np1) ,m,v,dumvar4,dumvar5]=sdot(m,v,1,v,1);    v(dumvar5~=v_orig)=dumvar5(dumvar5~=v_orig);
end;
wcnd = 1.;
nivn = fix(niv);
niv = 0;
%
if( y~=0.0 )
% **********************************************************************
for nr = 1 : n;
if( nivn==niv )
break;
end;
niv = fix(nr);
if( ix~=nr )
%
%     PIVOTING OF COLUMNS OF P MATRIX
%
nn = fix(n);
lix = fix(ix);
lr = fix(nr);
if( n~=ml18jr_17 )
nn = fix(ml18jr_17);
lix = fix(2.*ix - 1);
lr = fix(2.*nr - 1);
end;
if( nr~=1 )
kd = fix(lix - lr);
kj = fix(lr);
nrm1 = fix(lr - 1);
for j = 1 : nrm1;
psave = p(kj);
jk = fix(kj + kd);
p(kj) = p(jk);
p(jk) = psave;
kj = fix(kj + nn - j);
end; j = fix(nrm1+1);
jy = fix(jk + nmnr);
jz = fix(jy - kd);
p(jy) = p(jz);
end;
iz = fix(ip(lix));
ip(lix) = fix(ip(lr));
ip(lr) = fix(iz);
sv = s(ix);
s(ix) = s(nr);
s(nr) = sv;
if( n~=ml18jr_17 )
if( nr~=1 )
kj = fix(lr + 1);
for k = 1 : nrm1;
psave = p(kj);
jk = fix(kj + kd);
p(kj) = p(jk);
p(jk) = psave;
kj = fix(kj + ml18jr_17 - k);
end; k = fix(nrm1+1);
end;
iz = fix(ip(lix+1));
ip(lix+1) = fix(ip(lr+1));
ip(lr+1) = fix(iz);
end;
%
%     PIVOTING OF COLUMNS OF VECTORS
%
for l = 1 : m;
t = a(l,ix);
a(l,ix) = a(l,nr);
a(l,nr) = t;
end; l = fix(m+1);
end;
%
%     CALCULATE P(NR,NR) AS NORM SQUARED OF PIVOTAL VECTOR
%
jp = fix(jp + 1);
p(jp) = y;
ry = 1.0./y;
nmnr = fix(n - nr);
if( n~=ml18jr_17 )
nmnr = fix(ml18jr_17 -(2.*nr-1));
jp = fix(jp + 1);
p(jp) = 0.;
kp = fix(jp + nmnr);
p(kp) = y;
end;
if( nr~=n && nivn~=niv )
%
%    CALCULATE ORTHOGONAL PROJECTION VECTORS AND SEARCH FOR LARGEST NORM
%
y = 0.0;
ip1 = fix(nr + 1);
ix = fix(ip1);
%     ****************************************
for j = ip1 : n;
[dot ,m,dumvar2,dumvar4,dumvar4]=sdot(m,a(sub2ind(size(a),1,nr):end),1,a(sub2ind(size(a),1,j):end),1);      a(sub2ind(size(a),1,nr):end)=dumvar2; a(sub2ind(size(a),1,j):end)=dumvar4; 
jp = fix(jp + 1);
jq = fix(jp + nmnr);
if( n~=ml18jr_17 )
jq = fix(jq + nmnr - 1);
end;
p(jq) = p(jp) - dot.*(dot.*ry);
p(jp) = dot.*ry;
for i = 1 : m;
a(i,j) = a(i,j) - p(jp).*a(i,nr);
end; i = fix(m+1);
if( n~=ml18jr_17 )
kp = fix(jp + nmnr);
jp = fix(jp + 1);
pjp = ry.*prvec(m,a(sub2ind(size(a),1,nr):end),a(sub2ind(size(a),1,j):end));
p(jp) = pjp;
p(kp) = -pjp;
kp = fix(kp + 1);
p(kp) = ry.*dot;
for k = 1 : m2;
l = fix(m2 + k);
a(k,j) = a(k,j) - pjp.*a(l,nr);
a(l,j) = a(l,j) + pjp.*a(k,nr);
end; k = fix(m2+1);
p(jq) = p(jq) - pjp.*(pjp./ry);
end;
%
%     TEST FOR CANCELLATION IN RECURRENCE RELATION
%
if( p(jq)<=s(j).*ml5mco_2 )
[ p(jq) ,m,dumvar2,dumvar4,dumvar2]=sdot(m,a(sub2ind(size(a),1,j):end),1,a(sub2ind(size(a),1,j):end),1);      a(sub2ind(size(a),1,j):end)=dumvar2; 
end;
if( p(jq)>y )
y = p(jq);
ix = fix(j);
end;
end; j = fix(n+1);
if( n~=ml18jr_17 )
jp = fix(kp);
end;
%     ****************************************
if( ml18jr_11==1 )
ix = fix(ip1);
end;
%
%     RECOMPUTE NORM SQUARED OF PIVOTAL VECTOR WITH SCALAR PRODUCT
%
[y ,m,dumvar2,dumvar4,dumvar2]=sdot(m,a(sub2ind(size(a),1,ix):end),1,a(sub2ind(size(a),1,ix):end),1);      a(sub2ind(size(a),1,ix):end)=dumvar2; 
if( y<=ml5mco_3.*s(ix) )
iflag = 2;
wcnd = ml5mco_3;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
return;
end;
wcnd = min(wcnd,y./s(ix));
end;
%
%     COMPUTE ORTHOGONAL PROJECTION OF PARTICULAR SOLUTION
%
if( inhomo==1 )
lr = fix(nr);
if( n~=ml18jr_17 )
lr = fix(2.*nr - 1);
end;
w(lr) = sdot(m,a(sub2ind(size(a),1,nr):end),1,v,1).*ry;
for i = 1 : m;
v(i) = v(i) - w(lr).*a(i,nr);
end; i = fix(m+1);
if( n~=ml18jr_17 )
lr = fix(2.*nr);
w(lr) = ry.*prvec(m,v,a(sub2ind(size(a),1,nr):end));
for k = 1 : m2;
l = fix(m2 + k);
v(k) = v(k) + w(lr).*a(l,nr);
v(l) = v(l) - w(lr).*a(k,nr);
end; k = fix(m2+1);
end;
end;
end;
% **********************************************************************
%
%     TEST FOR LINEAR DEPENDENCE OF PARTICULAR SOLUTION
%
if( inhomo~=1 )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
return;
end;
if((n>1) &&(s(np1)<1.0) )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
return;
end;
v_orig=v;    [vnorm ,m,v,dumvar4,dumvar5]=sdot(m,v,1,v,1);    v(dumvar5~=v_orig)=dumvar5(dumvar5~=v_orig);
if( s(np1)~=0. )
wcnd = min(wcnd,vnorm./s(np1));
end;
if( vnorm>=ml5mco_3.*s(np1) )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
return;
end;
end;
iflag = 2;
wcnd = ml5mco_3;
else;
iflag = 1;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
return;
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
p_shape=zeros(p_shape);p_shape(:)=p(1:numel(p_shape));p=p_shape;
ip_shape=zeros(ip_shape);ip_shape(:)=ip(1:numel(ip_shape));ip=ip_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
end %subroutine mgsbv
%DECK MINFIT