(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



CG[A_,b_]:=Module[
{Tol=10^-8,n=Length[b],MaxIter=10^3,P,X,p,r,x,i=1},
P= ConstantArray[0,{MaxIter,n}];
X= ConstantArray[0,{MaxIter,n}];
(* Initializing the Conjugate Directions using b*)
r=b;
p=r/Sqrt[r.A.r];
x=(p.r)p;
{X[[i]],P[[i]]}={x,p};
While[
And[Norm[r]>Tol,i++<MaxIter],
(* Computing residual *)
r=b-A.x;
(* Conjugating *)
p=r;
p=p-(p.A.P[[i-1]]) P[[i-1]];
p=p/Sqrt[p.A.p];
(* Updating x *)
x+=(p.r)p;
(* Storing Data *)
{X[[i]],P[[i]]}={x,p}
];
i=Min[i,MaxIter];
{Take[X,i],Take[P,i]}
] 


BICG[A_,b_]:=Module[
{AT=Transpose[A],n=Length[b],
Tol=10^-8,MaxIter=10^3,
P,PHat,X,XHat,R, RHat,
p,pHat,r, rHat,x,k=1},
P=PHat=X= XHat=R=RHat= ConstantArray[0,{MaxIter,n}];
(* Initializing the Directions using Overscript[0, \[RightVector]] *)
p=pHat=r=rHat=b;
While[
And[Norm[r]>Tol,k++<MaxIter],
(* Computing residual *)
r=b-A.x;
(* Conjugating *)
p=r;
p=p-(p.A.P[[i-1]]) P[[i-1]];
p=p/Sqrt[p.A.p];
(* Updating x *)
x+=(p.r)p;
(* Storing Data *)
{X[[i]],P[[i]]}={x,p}
];
i=Min[i,MaxIter];
{Take[X,i],Take[P,i]}
] 



