[a,nrda,nrow,ncol,cols,colsav,rows,rowsav,anorm,scales,iscale,ic]=cscale(a,nrda,nrow,ncol,cols,colsav,rows,rowsav,anorm,scales,iscale,ic);

function [a,nrda,nrow,ncol,cols,colsav,rows,rowsav,anorm,scales,iscale,ic]=cscale(a,nrda,nrow,ncol,cols,colsav,rows,rowsav,anorm,scales,iscale,ic);
persistent alog2 ascale cs firstCall igo ip j k p s ten20 ten4 ; if isempty(firstCall),firstCall=1;end; 

if isempty(alog2), alog2=0; end;
if isempty(ascale), ascale=0; end;
if isempty(cs), cs=0; end;
if isempty(p), p=0; end;
if isempty(s), s=0; end;
if isempty(ten20), ten20=0; end;
if isempty(ten4), ten4=0; end;
if isempty(ip), ip=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(igo), igo=0; end;
%***BEGIN PROLOGUE  CSCALE
%***SUBSIDIARY
%***PURPOSE  Subsidiary to BVSUP
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (CSCALE-S, DCSCAL-D)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%     This routine scales the matrix A by columns when needed
%
%***SEE ALSO  BVSUP
%***ROUTINES CALLED  SDOT
%***REVISION HISTORY  (YYMMDD)
%   750601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%   910722  Updated AUTHOR section.  (ALS)
%***end PROLOGUE  CSCALE
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nrda])).*prod([nrda])-numel(a))],nrda,[]);
cols_shape=size(cols);cols=reshape(cols,1,[]);
colsav_shape=size(colsav);colsav=reshape(colsav,1,[]);
scales_shape=size(scales);scales=reshape(scales,1,[]);
rows_shape=size(rows);rows=reshape(rows,1,[]);
rowsav_shape=size(rowsav);rowsav=reshape(rowsav,1,[]);
%
if firstCall,   ten4 =[1.0e+4];  end;
if firstCall,  ten20=[1.0e+20];  end;
firstCall=0;
%
%***FIRST EXECUTABLE STATEMENT  CSCALE
igo=0;
if( iscale==(-1) )
%
if( ic~=0 )
for k = 1 : ncol;
[cols(k) ,nrow,dumvar2,dumvar4,dumvar2]=sdot(nrow,a(sub2ind(size(a),1,k):end),1,a(sub2ind(size(a),1,k):end),1);      a(sub2ind(size(a),1,k):end)=dumvar2; 
end; k = fix(ncol+1);
end;
%
ascale = anorm./ncol;
for k = 1 : ncol;
cs = cols(k);
if((cs>ten4.*ascale) ||(ten4.*cs<ascale) )
igo=1;
break;
end;
if((cs<1../ten20) ||(cs>ten20) )
igo=1;
break;
end;
end;
end;
%
if(igo==0)
for k = 1 : ncol;
scales(k) = 1.;
end; k = fix(ncol+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
cols_shape=zeros(cols_shape);cols_shape(:)=cols(1:numel(cols_shape));cols=cols_shape;
colsav_shape=zeros(colsav_shape);colsav_shape(:)=colsav(1:numel(colsav_shape));colsav=colsav_shape;
scales_shape=zeros(scales_shape);scales_shape(:)=scales(1:numel(scales_shape));scales=scales_shape;
rows_shape=zeros(rows_shape);rows_shape(:)=rows(1:numel(rows_shape));rows=rows_shape;
rowsav_shape=zeros(rowsav_shape);rowsav_shape(:)=rowsav(1:numel(rowsav_shape));rowsav=rowsav_shape;
return;
end;
%
alog2 = log(2.);
anorm = 0.;
for k = 1 : ncol;
cs = cols(k);
if( cs~=0. )
p = log(cs)./alog2;
ip = fix(-0.5.*p);
s = 2..^ip;
scales(k) = s;
if( ic~=1 )
cols(k) = s.*s.*cols(k);
anorm = anorm + cols(k);
colsav(k) = cols(k);
end;
for j = 1 : nrow;
a(j,k) = s.*a(j,k);
end; j = fix(nrow+1);
else;
scales(k) = 1.;
end;
end; k = fix(ncol+1);
%
if( ic==0 )
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
cols_shape=zeros(cols_shape);cols_shape(:)=cols(1:numel(cols_shape));cols=cols_shape;
colsav_shape=zeros(colsav_shape);colsav_shape(:)=colsav(1:numel(colsav_shape));colsav=colsav_shape;
scales_shape=zeros(scales_shape);scales_shape(:)=scales(1:numel(scales_shape));scales=scales_shape;
rows_shape=zeros(rows_shape);rows_shape(:)=rows(1:numel(rows_shape));rows=rows_shape;
rowsav_shape=zeros(rowsav_shape);rowsav_shape(:)=rowsav(1:numel(rowsav_shape));rowsav=rowsav_shape;
return;
end;
%
for k = 1 : nrow;
nrda_orig=nrda;    [rows(k) ,ncol,dumvar2,nrda,dumvar2,dumvar6]=sdot(ncol,a(sub2ind(size(a),k,1):end),nrda,a(sub2ind(size(a),k,1):end),nrda);    nrda(dumvar6~=nrda_orig)=dumvar6(dumvar6~=nrda_orig);      a(sub2ind(size(a),k,1):end)=dumvar2; 
rowsav(k) = rows(k);
anorm = anorm + rows(k);
end; k = fix(nrow+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
cols_shape=zeros(cols_shape);cols_shape(:)=cols(1:numel(cols_shape));cols=cols_shape;
colsav_shape=zeros(colsav_shape);colsav_shape(:)=colsav(1:numel(colsav_shape));colsav=colsav_shape;
scales_shape=zeros(scales_shape);scales_shape(:)=scales(1:numel(scales_shape));scales=scales_shape;
rows_shape=zeros(rows_shape);rows_shape(:)=rows(1:numel(rows_shape));rows=rows_shape;
rowsav_shape=zeros(rowsav_shape);rowsav_shape(:)=rowsav(1:numel(rowsav_shape));rowsav=rowsav_shape;
end
%DECK CSCAL