template<class T>
void BandMtx<T>::GaussElimPP(Vcr<T>&bb) const{
	if (nrows != bb.size())
		error("matrix-vector sizes not match");
	BandMtx<T> tx(nrows,bwleft,min(nrows-1,bwleft+bwrit));
	 for(int i=0;i<nrows;i++)
       for(int j=-bwleft;j<=bwrit;j++)
		   tx[i][j]=bdmx[i][j];
	   Vcr<int> pvt(nrows);  //保存主元信息
	     
	   //部分列主元的LU分解
	   const int nrowsmone=tx.nrow-1;
	   for (int k=0;k<nrowsmone-k;k++){
        int kbrow=min(nrowsmone-k,tx.bwleft);
		int kbcol=min(nrowsmone-k,tx.bwrit);

		//找第k列的主元
		int pc=k;
		double aet=abs(tx[k][0]);
		for (int i=1;i<=kbrow;i++)
		{if (abs(tx[k+i][-i])>aet)
		{ aet=abs(tx[k+i][-i]);
		    pc=k+i;
		}
		}
		if(!aet) error("pivot is zero in GaussElimPP");
		pvt[k]=pc;


		//交換U,而不是L的主元行和第k行
		if(pc!=k){
			for(int j=0;j<=kbcol;j++)
				swap(tx[pc][k+j-pc],tx[k][j]);
		}


		//接下來消去該列元素
		for(int i=1;i<=kbrow;i++){
			int kpi=k+i;
			if (tx[kpi][-i]!=0){
				T dmul=tx[kpi][-i]/tx[k][0];
                tx[kpi][-i]=dmul;
				for(int j=1;j<=kbcol;j++)
					tx[kpi][j-i]-=dmul*tx[k][j];
			}
		}
}
  


//向前代入法求解LY=Pb
for(int k=0; k<nrowsmone;k++){
	int kbrow=min(nrowsmone-k,tx.bwleft);
	int pvtk=pvt[k];
	T sb=bb[pvtk];
	if(k!=pvtk)
		swap(bb[k],bb[pvtk]);
	 for(int j=1;j<=kbrow;j++)
		 bb[k+j]-=tx[k+j][-j]*sb;
}


//向后代入法求解Ux=y
for(int k=nrowsmone;k>=0;k--){
	int kb=min(nrowsmone-k,tx.bwrit);
	for(int j=1;j<=kb;j++) bb[k]-=tx[k][j]*bb[k+j];
	bb[k]/=tx[k][0];
}
}   //GaussElimpp()結束