function [ x, istop, itn, nrmgrad] = ...         truncLS3(VEC, m, n, aprodname, mu, D, b, damp, atol, btol, conlim, ...	 itnlim,w1,b1,a1,deriv1,df1,w2,b2,a2,deriv2,df2,w3,b3,a3,deriv3,df3,P,T)%        [ x, istop, itn, nrmgrad ] = ...%        truncLS3(VEC, m, n, aprodname, mu, D, b, damp, atol, btol, conlim, itnlim,...%	 w1,b1,a1,deriv1,w2,b2,a2,deriv2,w3,b3,a3,deriv3,P,T)%%truncLS3 uses a truncation technique to approximately solve% (   A)D*D(inv)x = b   in the least squares sence if damp = 0 and m>=n.% (mu*I)% where D is diagonal scaling matrix.% D holds the diagonal elements.% 		Ax = b  if m = n and damp = 0,% or           min || (b)  -  (  A   )x ||   otherwise.%                  || (0)     (damp I)  ||2% A  is an m by n matrix defined by  y = aprod( mode,m,n,x,mu,D),% where the parameter 'aprodname' refers to a function 'aprod' that% performs the matrix-vector operations.% If mode = 1,   aprod  must return  y = Ax   without altering x.% If mode = 2,   aprod  must return  y = A'x  without altering x.% WARNING:   The file containing the function 'aprod'%            must not be called aprodname.m !!!!% LSQR uses an iterative (conjugate-gradient-like) method.% For further information, see C. C. Paige and M. A. Saunders:% 1. LSQR: An algorithm for sparse linear equations and sparse least squares,%    ACM TOMS 8, 1 (Mar 1982) 43-71.% 2. Algorithm 583.  LSQR: Sparse linear equations and least squares problems,%    ACM TOMS 8, 2 (Jun 1982) 195-209.%% Other potential output parameters:%    anorm, acond, rnorm, arnorm, xnorm%%        1990: Derived from Fortran 77 version of LSQR.% 22 May 1992: bbnorm was used incorrectly.  Replaced by anorm.% 26 Oct 1992: More input and output parameters added.% 01 Sep 1994: Matrix-vector routine is now a parameter 'aprodname'.%              Print log reformatted.%     Initialize.msg=['The exact solution is  x = 0                              '     'Ax - b is small enough, given atol, btol                  '     'The least-squares solution is good enough, given atol     '     'The estimate of cond(Abar) has exceeded conlim            '     'Ax - b is small enough for this machine                   '     'The least-squares solution is good enough for this machine'     'Cond(Abar) seems to be too large for this machine         '     'The iteration limit has been reached                      '];itn    = 0;		istop  = 0;		nstop  = 0;ctol   = 0;		if conlim > 0, ctol = 1/conlim; end;anorm  = 0;		acond  = 0;dampsq = damp^2;	ddnorm = 0;		res2   = 0;xnorm  = 0;		xxnorm = 0;		z      = 0;cs2    = -1;		sn2    = 0;major_it=VEC(1);termind=VEC(2);% Set up the first vectors u and v for the bidiagonalization.% These satisfy  beta*u = b,  alfa*v = A'u.u      = b(1:m);	x    = zeros(n,1);alfa   = 0;		beta = norm( u );if beta > 0   u = (1/beta) * u;	v = D .* feval( aprodname, 2, m, n, u, mu, D, ...   w1,b1,a1,deriv1,df1,w2,b2,a2,deriv2,df2,w3,b3,a3,deriv3,df3,P,T);   alfa = norm( v );%**************nrmgrad=alfa*beta;eta=min(1/major_it, nrmgrad);%nrmrhs=alfa;endif alfa > 0   v = (1/alfa) * v;    w = v;endarnorm = alfa * beta;   if arnorm == 0, disp(msg(1,:)); return, endrhobar = alfa;		phibar = beta;		bnorm  = beta;rnorm  = beta;head1  = '   Itn      x(1)      Function';head2  = ' Compatible   LS      Norm A   Cond A';%disp([head1 head2])test1  = 1;		test2  = alfa / beta;str1   = sprintf( '%6g %12.5e %10.3e',   itn, x(1), rnorm );str2   = sprintf( '  %8.1e %8.1e',       test1, test2 );%disp([str1 str2])%------------------------------------------------------------------%     Main iteration loop.%------------------------------------------------------------------while itn < itnlim      itn = itn + 1;%     Perform the next step of the bidiagonalization to obtain the%     next  beta, u, alfa, v.  These satisfy the relations%                beta*u  =  a*v   -  alfa*u,%                alfa*v  =  A'*u  -  beta*v.%***************      u    = feval( aprodname, 1, m, n, D .* v, mu, D, ...      w1,b1,a1,deriv1,df1,w2,b2,a2,deriv2,df2,w3,b3,a3,deriv3,df3,P,T)  -  alfa*u;      beta = norm( u );      if beta > 0         u     = (1/beta) * u;         anorm = norm([anorm alfa beta damp]);%****************         v     = D .* feval( aprodname, 2, m, n, u, mu, D, ...	 w1,b1,a1,deriv1,df1,w2,b2,a2,deriv2,df2,w3,b3,a3,deriv3,df3,P,T)  -  beta*v;         alfa  = norm( v );         if alfa > 0,  v = (1/alfa) * v; end      end%     Use a plane rotation to eliminate the damping parameter.%     This alters the diagonal (rhobar) of the lower-bidiagonal matrix.      rhobar1 = norm([rhobar damp]);      cs1     = rhobar / rhobar1;      sn1     = damp   / rhobar1;      psi     = sn1 * phibar;      phibar  = cs1 * phibar;%     Use a plane rotation to eliminate the subdiagonal element (beta)%     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.      rho     = norm([rhobar1 beta]);      cs      =   rhobar1/ rho;      sn      =   beta   / rho;      theta   =   sn * alfa;      rhobar  = - cs * alfa;      phi     =   cs * phibar;      phibar  =   sn * phibar;      tau     =   sn * phi;%     Update x and w.      ddnorm  = ddnorm + (( 1 / rho)*norm( w ))^2;      x       = x      + (phi / rho)*w;      w       = v      - (theta/rho)*w;%     Use a plane rotation on the right to eliminate the%     super-diagonal element (theta) of the upper-bidiagonal matrix.%     Then use the result to estimate  norm(x).      delta   =   sn2 * rho;      gambar  = - cs2 * rho;      rhs     =   phi  -  delta * z;      zbar    =   rhs / gambar;      xnorm   =   sqrt(xxnorm + zbar^2);      gamma   =   norm([gambar theta]);      cs2     =   gambar / gamma;      sn2     =   theta  / gamma;      z       =   rhs    / gamma;      xxnorm  =   xxnorm  +  z^2;%     Test for convergence.%     First, estimate the condition of the matrix  Abar,%     and the norms of  rbar  and  Abar'rbar.      acond   =   anorm * sqrt( ddnorm );      res1    =   phibar^2;      res2    =   res2  +  psi^2;      rnorm   =   sqrt( res1 + res2 );      arnorm  =   alfa * abs( tau );%     Now use these norms to estimate certain other quantities,%     some of which will be small near a solution.      test1   =   rnorm / bnorm;RNORM(itn,1)=rnorm;RNORM(itn,2)=test1;RNORM(itn,3)=xnorm;RNORM(itn,4)=arnorm;RNORM(itn,5)=arnorm/bnorm;RNORM(itn,6)=arnorm/nrmgrad;RNORM(itn,7)=eta;      test2   =   arnorm/( anorm * rnorm );      test3   =       1 / acond;      t1      =   test1 / (1    +  anorm * xnorm / bnorm);      rtol    =   btol  +  atol *  anorm * xnorm / bnorm;%     The following tests guard against extremely small values of%     atol, btol  or  ctol.  (The user may have set any or all of%     the parameters  atol, btol, conlim  to 0.)%     The effect is equivalent to the normal tests using%     atol = eps,  btol = eps,  conlim = 1/eps.%if (RNORM(itn,5) < 1e-3) & (RNORM(itn,3) < 1e-2),   istop = 8; end%if RNORM(itn,5) < 0.5, istop = 8; endif termind == 1	if arnorm < nrmgrad*eta, istop = 9; end	if itn >= itnlim,   istop = 7; endelseif termind == 3	if itn >= itnlim,   istop = 7; endelse %termind == 2	if itn >= itnlim,   istop = 7; end      if 1 + test3  <= 1, istop = 6; end      if 1 + test2  <= 1, istop = 5; end      if 1 + t1     <= 1, istop = 4; end%     Allow for tolerances set by the user.      if  test3 <= ctol,  istop = 3; end      if  test2 <= atol,  istop = 2; end      if  test1 <= rtol,  istop = 1; endendif istop > 0, break, end%     See if it is time to print something.      prnt = 0;      if n     <= 40       , prnt = 1; end      if itn   <= 10       , prnt = 1; end      if itn   >= itnlim-10, prnt = 1; end      if rem(itn,10) == 0  , prnt = 1; end      if test3 <=  2*ctol  , prnt = 1; end      if test2 <= 10*atol  , prnt = 1; end      if test1 <= 10*rtol  , prnt = 1; end      if istop ~=  0       , prnt = 1; end      if prnt == 1         str1 = sprintf( '%6g %12.5e %10.3e',   itn, x(1), rnorm );         str2 = sprintf( '  %8.1e %8.1e',       test1, test2 );         str3 = sprintf( ' %8.1e %8.1e',        anorm, acond );         %disp([str1 str2 str3])      end      if istop > 0, break, endendx=D .* x;%[istop itn bnorm nrmgrad]%RNORM%pause%     End of iteration loop.%     Print the stopping condition.%disp(' ')%disp('LSQR finished')%disp(msg(istop+1,:))%disp(' ')str1 = sprintf( 'istop  =%8g   itn    =%8g',      istop, itn    );str2 = sprintf( 'anorm  =%8.1e   acond  =%8.1e',  anorm, acond  );str3 = sprintf( 'rnorm  =%8.1e   arnorm =%8.1e',  rnorm, arnorm );str4 = sprintf( 'bnorm  =%8.1e   xnorm  =%8.1e',  bnorm, xnorm  );%disp([str1 '   ' str2])%disp([str3 '   ' str4])%     End of LSQR