/*
description:
	fuer genaue aufgabenstellung siehe aufgabenblatt 2
tab stops:
	2 (otherwise hardly readable!)
*/

#include "matrix.h"
#include <math.h>
#include <float.h>
#include <stdio.h>
#include <string.h>

void aufgabe6(vektor&, matrix&, double=0, int=0, double=0);
int ivektoriteration(vektor&, matrix&, int=1, double=0, int=0, double=0);

int ivektoriteration(vektor &start, matrix &A, int verbose, double ew, int eigenwert_flag, double lambda){
	matrix L(A.get_dim_zeilen(), A.get_dim_spalten());
	matrix R(A.get_dim_zeilen(), A.get_dim_spalten());
	if((A-lambda).LR(L, R)==0) cout <<"error: matrix A-lambda nicht LR-zerlegbar. ab jetzt kommt bloedsinn heraus...";
	if(verbose){
		cout <<"  L="<<endl;
		L.ausgabe();
		cout <<"  R="<<endl;
		R.ausgabe();
		if(eigenwert_flag) cout <<endl<<"  eigenwert (bekannt) = "<<ew<<endl;
		cout <<"  startvektor="<<start<<endl;
		cout.unsetf(ios::scientific);	// scientific-notation abschalten
		cout <<"  iteration:"<<endl;
		cout <<"  schritte  |mu^(k)|   E_rel^(k)   kappa^(k)"<<endl;
	}
	int k=0;
	int dim=start.get_dim();
	double mu=0, mu_old, E_rel, kappa, d=0, d_k1=0, d_k2;
	vektor u(start), v(start), w(dim);
	const double TOL=1.0e-4;
	while(k<2 || (d>=TOL && k<1000)){
		if(verbose) cout <<"     "<<((k<10)?" ":"")<<k;
		w=L.solve_linksunten(u);
		v=R.solve_rechtsoben(w);
		if(v.norm()>0) u=v*(1.0/v.norm());
		else {
			cout <<"error: division by zero, breche schleife nach "<<k<<" it-schritten ab"<<endl;
			break;
		}
		if(verbose){
			cout <<"      ";
			printf("%5.4f", fabs(mu));
			if(eigenwert_flag){
				E_rel=fabs(mu-ew)/fabs(ew);
				cout <<"     ";
				printf("%5.4f", E_rel);
			}
		}
		mu_old=mu;
		mu=v.norm();
		d_k2=d_k1;
		d_k1=d;
		d=fabs(mu-mu_old);
		if(verbose && k>1){
			if(d_k2==0 || d_k1==0 || d_k1==d_k2) cout <<"error: division by zero (berechnung von kappa)"<<endl;
			else{
				kappa=log(d/d_k1)/log(d_k1/d_k2);
				cout <<"      ";
				printf("%5.4f", kappa);
			}
		}
		if(verbose) cout <<endl;
		k++;
	}
	if(verbose){
		cout <<endl<<"  abbruch nach "<<--k<<" schritten."<<endl;
		cout <<"  iterierte werte:"<<endl;
		cout <<"  lambda_s="<<(1/mu)+lambda<<endl;
		cout <<"  zug. EV s_1="<<u;
	}
	return k;
}

void aufgabe6(vektor &start, matrix &A, double ew, int eigenwert_flag, double lambda){
	cout <<"6.aufg:"<<endl;
	cout <<"======="<<endl;
	ivektoriteration(start, A, 1, ew, eigenwert_flag, lambda);
}

int main(){
	cout <<endl;
	int dim=10;
	double lambda=0.65;// shift
	double lambda_s=6.9028e-1;// bekannter EW
	matrix A(dim, dim); // matrix
	A.tridiag(-1, 2, -1);
	vektor u(dim); // startvektor
	u.set(0, 1);
	aufgabe6(u, A, lambda_s, 1, lambda);
	cout <<endl;
	return 0;
}