[dqdotiresult,n,db,qc,dx,incx,dy,incy]=dqdoti(n,db,qc,dx,incx,dy,incy);

function [dqdotiresult,n,db,qc,dx,incx,dy,incy]=dqdoti(n,db,qc,dx,incx,dy,incy);
dqdotiresult=[];
persistent firstCall i i1 ix iy qx qy ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(i1), i1=0; end;
if isempty(ix), ix=0; end;
if isempty(iy), iy=0; end;
global mpcom_1; if isempty(mpcom_1), mpcom_1=0; end;
global mpcom_4; if isempty(mpcom_4), mpcom_4=0; end;
global mpcom_3; if isempty(mpcom_3), mpcom_3=0; end;
global mpcom_5; if isempty(mpcom_5), mpcom_5=0; end;
global mpcom_6; if isempty(mpcom_6), mpcom_6=zeros(1,30); end;
global mpcom_2; if isempty(mpcom_2), mpcom_2=0; end;
%***BEGIN PROLOGUE  DQDOTI
%***PURPOSE  Compute the inner product of two vectors with extended
%            precision accumulation and result.
%***LIBRARY   SLATEC
%***CATEGORY  D1A4
%***TYPE      doubleprecision (DQDOTI-D)
%***KEYWORDS  DOT PRODUCT, INNER PRODUCT
%***AUTHOR  Lawson, C. L., (JPL)
%           Hanson, R. J., (SNLA)
%           Kincaid, D. R., (U. of Texas)
%           Krogh, F. T., (JPL)
%***DESCRIPTION
%
%                B L A S  Subprogram
%    Description of parameters
%
%     --Input--
%        N  number of elements in input vector(s)
%       DB  doubleprecision scalar to be added to inner product
%       QC  extended precision scalar to be added
%       DX  doubleprecision vector with N elements
%     INCX  storage spacing between elements of DX
%       DY  doubleprecision vector with N elements
%     INCY  storage spacing between elements of DY
%
%     --Output--
%   DQDOTI  doubleprecision result
%       QC  extended precision result
%
%     D.P. dot product with extended precision accumulation (and result)
%     QC and DQDOTI are set = DB + sum for I = 0 to N-1 of
%       DX(LX+I*INCX) * DY(LY+I*INCY),  where QC is an extended
%       precision result which can be used as input to DQDOTA,
%       and LX = 1 if INCX .GE. 0, else LX = (-INCX)*N, and LY is
%       defined in a similar way using INCY.  The MP package by
%       Richard P. Brent is used for the extended precision arithmetic.
%
%     Fred T. Krogh,  JPL,  1977,  June 1
%
%     The common block for the MP package is named MPCOM.  If local
%     variable I1 is zero, DQDOTI calls MPBLAS to initialize the MP
%     package and reset I1 to 1.
%
%    The argument QC(*), and the local variables QX and QY are INTEGER
%    arrays of size 30.  See the comments in the routine MPBLAS for the
%    reason for this choice.
%
%***REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
%                 Krogh, Basic linear algebra subprograms for Fortran
%                 usage, Algorithm No. 539, Transactions on Mathematical
%                 Software 5, 3 (September 1979), pp. 308-323.
%***ROUTINES CALLED  MPADD, MPBLAS, MPCDM, MPCMD, MPMUL
%***COMMON BLOCKS    MPCOM
%***REVISION HISTORY  (YYMMDD)
%   791001  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%   930124  Increased Array sizes for SUN -r8.  (RWC)
%***end PROLOGUE  DQDOTI
dx_shape=size(dx);dx=reshape(dx,1,[]);
dy_shape=size(dy);dy=reshape(dy,1,[]);
if isempty(qx), qx=zeros(1,30); end;
if isempty(qy), qy=zeros(1,30); end;
% common :: ;
%% common /mpcom / mpb , mpt , mpm , mplun , mpmxr , mpr(30);
%% common /mpcom / mpcom_1 , mpcom_2 , mpcom_3 , mpcom_4 , mpcom_5 , mpcom_6(30);
if firstCall,   i1=[0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DQDOTI
if( i1==0 )
[i1]=mpblas(i1);
end;
qc(1) = 0;
if( db~=0.0d0 )
[db,qx]=mpcdm(db,qx);
qc_orig=qc;    [qc,qx,dumvar3]=mpadd(qc,qx,qc);    qc(dumvar3~=qc_orig)=dumvar3(dumvar3~=qc_orig);
end;
if( n~=0 )
ix = 1;
iy = 1;
if( incx<0 )
ix =fix((-n+1).*incx + 1);
end;
if( incy<0 )
iy =fix((-n+1).*incy + 1);
end;
for i = 1 : n;
[dx(ix),qx]=mpcdm(dx(ix),qx);
[dy(iy),qy]=mpcdm(dy(iy),qy);
qx_orig=qx;    [qx,qy,dumvar3]=mpmul(qx,qy,qx);    qx(dumvar3~=qx_orig)=dumvar3(dumvar3~=qx_orig);
qc_orig=qc;    [qc,qx,dumvar3]=mpadd(qc,qx,qc);    qc(dumvar3~=qc_orig)=dumvar3(dumvar3~=qc_orig);
ix = fix(ix + incx);
iy = fix(iy + incy);
end; i = fix(n+1);
end;
[qc,dqdotiresult]=mpcmd(qc,dqdotiresult);
dx_shape=zeros(dx_shape);dx_shape(:)=dx(1:numel(dx_shape));dx=dx_shape;
dy_shape=zeros(dy_shape);dy_shape(:)=dy(1:numel(dy_shape));dy=dy_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',qc); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',n); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',incy); evalin('caller',[inputname(7),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',incx); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(6)), assignin('caller','FUntemp',dy); evalin('caller',[inputname(6),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',dx); evalin('caller',[inputname(4),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',db); evalin('caller',[inputname(2),'=FUntemp;']); end
end
%DECK DQELG