| [transa,transb,m,n,kk,alpha,a,lda,b,ldb,beta,c,ldc,ct,g,cc,ldcc,eps,err,ftl,nout,mv,kprint]=cmmch(transa,transb,m,n,kk,alpha,a,lda,b,ldb,beta,c,ldc,ct,g,cc,ldcc,eps,err,ftl,nout,mv,kprint); |
function [transa,transb,m,n,kk,alpha,a,lda,b,ldb,beta,c,ldc,ct,g,cc,ldcc,eps,err,ftl,nout,mv,kprint]=cmmch(transa,transb,m,n,kk,alpha,a,lda,b,ldb,beta,c,ldc,ct,g,cc,ldcc,eps,err,ftl,nout,mv,kprint);
persistent cl ctrana ctranb erri i j k rone rzero trana tranb zero ;
if isempty(zero), zero=complex(0.0,0.0) ; end;
if isempty(rzero), rzero=0.0; end;
if isempty(rone), rone=1.0 ; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
b_shape=size(b);b=reshape([b(:).',zeros(1,ceil(numel(b)./prod([ldb])).*prod([ldb])-numel(b))],ldb,[]);
c_shape=size(c);c=reshape([c(:).',zeros(1,ceil(numel(c)./prod([ldc])).*prod([ldc])-numel(c))],ldc,[]);
cc_shape=size(cc);cc=reshape([cc(:).',zeros(1,ceil(numel(cc)./prod([ldcc])).*prod([ldcc])-numel(cc))],ldcc,[]);
ct_shape=size(ct);ct=reshape(ct,1,[]);
g_shape=size(g);g=reshape(g,1,[]);
if isempty(cl), cl=0; end;
if isempty(erri), erri=0; end;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(ctrana), ctrana=false; end;
if isempty(ctranb), ctranb=false; end;
if isempty(trana), trana=false; end;
if isempty(tranb), tranb=false; end;
% intrinsic abs , aimag , conjg , max , real :: ;
% abs1= @(cl) abs(real(cl)) + abs(aimag(cl));real :: abs1;
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
trana = strcmp(deblank(transa),deblank('T')) | strcmp(deblank(transa),deblank('C'));
tranb = strcmp(deblank(transb),deblank('T')) | strcmp(deblank(transb),deblank('C'));
ctrana = strcmp(deblank(transa),deblank('C'));
ctranb = strcmp(deblank(transb),deblank('C'));
for j = 1 : n;
for i = 1 : m;
ct(i) = zero;
g(i) = rzero;
end; i = fix(m+1);
if( ~trana && ~tranb )
for k = 1 : kk;
for i = 1 : m;
ct(i) = ct(i) + a(i,k).*b(k,j);
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = g(i) + abs1(a(i,k)).*abs1(b(k,j));
end; i = fix(m+1);
end; k = fix(kk+1);
elseif( trana && ~tranb ) ;
if( ctrana )
for k = 1 : kk;
for i = 1 : m;
ct(i) = ct(i) + conj(a(k,i)).*b(k,j);
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = g(i) + abs1(a(k,i)).*abs1(b(k,j));
end; i = fix(m+1);
end; k = fix(kk+1);
else;
for k = 1 : kk;
for i = 1 : m;
ct(i) = ct(i) + a(k,i).*b(k,j);
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = g(i) + abs1(a(k,i)).*abs1(b(k,j));
end; i = fix(m+1);
end; k = fix(kk+1);
end;
elseif( ~trana && tranb ) ;
if( ctranb )
for k = 1 : kk;
for i = 1 : m;
ct(i) = ct(i) + a(i,k).*conj(b(j,k));
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = g(i) + abs1(a(i,k)).*abs1(b(j,k));
end; i = fix(m+1);
end; k = fix(kk+1);
else;
for k = 1 : kk;
for i = 1 : m;
ct(i) = ct(i) + a(i,k).*b(j,k);
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = g(i) + abs1(a(i,k)).*abs1(b(j,k));
end; i = fix(m+1);
end; k = fix(kk+1);
end;
elseif( trana && tranb ) ;
if( ctrana )
if( ctranb )
for k = 1 : kk;
for i = 1 : m;
ct(i) = ct(i) + conj(a(k,i)).*conj(b(j,k));
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = g(i) + abs1(a(k,i)).*abs1(b(j,k));
end; i = fix(m+1);
end; k = fix(kk+1);
else;
for k = 1 : kk;
for i = 1 : m;
ct(i) = ct(i) + conj(a(k,i)).*b(j,k);
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = g(i) + abs1(a(k,i)).*abs1(b(j,k));
end; i = fix(m+1);
end; k = fix(kk+1);
end;
elseif( ctranb ) ;
for k = 1 : kk;
for i = 1 : m;
ct(i) = ct(i) + a(k,i).*conj(b(j,k));
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = g(i) + abs1(a(k,i)).*abs1(b(j,k));
end; i = fix(m+1);
end; k = fix(kk+1);
else;
for k = 1 : kk;
for i = 1 : m;
ct(i) = ct(i) + a(k,i).*b(j,k);
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = g(i) + abs1(a(k,i)).*abs1(b(j,k));
end; i = fix(m+1);
end; k = fix(kk+1);
end;
end;
for i = 1 : m;
ct(i) = alpha.*ct(i) + beta.*c(i,j);
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
g(i) = abs1(alpha).*g(i) + abs1(beta).*abs1(c(i,j));
end; i = fix(m+1);
err = zero;
for i = 1 : m;
abs1= @(cl) abs(real(cl)) + abs(imag(cl));
erri = abs1(ct(i)-cc(i,j))./eps;
if( g(i)~=rzero )
erri = erri./g(i);
end;
err = max(err,erri);
if( err.*sqrt(eps)>=rone )
ftl = true;
if( kprint>=2 )
writef(nout,[' ** FATAL ERROR - COMPUTED RESULT IS LESS THAN HAL','F ACCURATE **', '\n ' ,' EXPECTED RE','SULT COMPUTED RESULT' ' \n']);
for k = 1 : m;
if( mv )
writef(nout,[repmat(' ',1,1),'%7i',repmat([' (','%15.6f',',','%15.6f',')'] ,1,2) ' \n'], k , ct(k) , cc(k,j));
else;
writef(nout,[repmat(' ',1,1),'%7i',repmat([' (','%15.6f',',','%15.6f',')'] ,1,2) ' \n'], k , cc(k,j) , ct(k));
end;
end; k = fix(m+1);
end;
end;
end; i = fix(m+1);
end; j = fix(n+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
cc_shape=zeros(cc_shape);cc_shape(:)=cc(1:numel(cc_shape));cc=cc_shape;
ct_shape=zeros(ct_shape);ct_shape(:)=ct(1:numel(ct_shape));ct=ct_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
return;
%format (' ** FATAL ERROR - COMPUTED RESULT IS LESS THAN HAL','F ACCURATE **',/' EXPECTED RE','SULT COMPUTED RESULT');
%format (1X,i7,2(' (',g15.6,',',g15.6,')'));
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
cc_shape=zeros(cc_shape);cc_shape(:)=cc(1:numel(cc_shape));cc=cc_shape;
ct_shape=zeros(ct_shape);ct_shape(:)=ct(1:numel(ct_shape));ct=ct_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
end %subroutine cmmch
|
|