function [Af] =aprod3aug(mode,m,n,f,mu,rw, ...	  w1,b1,Output_Layer2,Deriv_Layer2,df1,w2,b2,Output_Layer3,...	  Deriv_Layer3,df2,w3,b3,Output_Layer4,Deriv_Layer4,df3,P,T)%Call: [Af] =aprod3aug(mode,m,n,f,mu,rw, ...%	  w1,b1,Output_Layer2,Deriv_Layer2,df1,w2,b2,Output_Layer3,Deriv_Layer3,%	  ..., df2,w3,b3,Output_Layer4,Deriv_Layer4,df3,P,T)%This routine computes Jaug*f or Jaug'*f, where J is the Jacobian corresponding%to a FFANN with two hidden layers.%Jaug = ( J  )%       (mu*I)%It should normally be used with the routine truncLS3 which is a%modification of LSQR written by Paige and Saunders%% If mode == 1 aprod3aug returns y=Jaug*f% If mode == 2 aprod3aug returns y=Jaug'*f%[Fst,q]=size(P);	% # nodes in Input layer (first layer)[Snd,r]=size(w1);	% # nodes in 1'st hidden layer (second layer)if Fst ~= r		% q is # samples presented to the inputs   error('In aprod3aug: Fst .ne. r');end[Thrd,r]=size(w2);	% # nodes in 2'nd hidden layer (third layer)if Snd ~= r   error('In aprod3aug: Snd .ne. r')end[Fth,r]=size(w3);	% # nodes in output layer (forth layer)Om=Fst*Snd+Snd;		% # weights in first hidden layerNm=Om+Snd*Thrd+Thrd;	% # weights up to and including second hidden layerMm=Nm+Thrd*Fth+Fth;	% # total weightslenf=max(size(f));if mode == 1  if Mm ~= lenf     error('In aprod3aug: Mm .ne. lenf (mode == 1)');  else     Af=zeros(q*Fth+Mm,1);  endelseif (q*Fth+Mm) ~= lenf      error('In aprod3aug: q*Fth+Mm .ne. lenf (mode == 2)');    else      Af=zeros(Mm,1);end%if mode == 1%Compute Jaug*f by using forward mode%Unpack the vector elements such that same kind of computation as%in the function evaluation can be done.%First hidden layerTEMP=ones(1,q);w1=reshape(f(1:Fst*Snd),Snd,Fst);b1=f(Fst*Snd+1:Fst*Snd+Snd);O_L2=b1*TEMP+w1*P;S2=Deriv_Layer2 .* O_L2;%Second hidden layerw2temp=reshape(f(Om+1:Om+Snd*Thrd),Thrd,Snd);b2=f(Om+Snd*Thrd+1:Om+Snd*Thrd+Thrd);O_L3=b2*TEMP+w2*S2+w2temp*Output_Layer2;S3=Deriv_Layer3 .* O_L3;%Output layerw3temp=reshape(f(Nm+1:Nm+Fth*Thrd),Fth,Thrd);b3=f(Nm+Fth*Thrd+1:Mm);O_L4=b3*TEMP+w3*S3+w3temp*Output_Layer3;S4=Deriv_Layer4 .* O_L4;Af=[S4(:); mu*f];else % Now mode should be 2   %Compute Jaug'*f   e=reshape(f(1:Fth*q),Fth,q);   d3=feval(df3,Output_Layer4,e);   d2=feval(df2,Output_Layer3,d3,w3);   d1=feval(df1,Output_Layer2,d2,w2);   [dw1,db1]=learnbp(P,d1,1);   [dw2,db2]=learnbp(Output_Layer2,d2,1);   [dw3,db3]=learnbp(Output_Layer3,d3,1);   Af=[dw1(:); db1(:); dw2(:); db2(:); dw3(:); db3(:)];   Af=Af+mu*f(q*Fth+1:q*Fth+Mm);end %of if mode == 1end