/*******************************************************************************
* Copyright 2005-2016 Intel Corporation All Rights Reserved.
*
* The source code,  information  and material  ("Material") contained  herein is
* owned by Intel Corporation or its  suppliers or licensors,  and  title to such
* Material remains with Intel  Corporation or its  suppliers or  licensors.  The
* Material  contains  proprietary  information  of  Intel or  its suppliers  and
* licensors.  The Material is protected by  worldwide copyright  laws and treaty
* provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
* modified, published,  uploaded, posted, transmitted,  distributed or disclosed
* in any way without Intel's prior express written permission.  No license under
* any patent,  copyright or other  intellectual property rights  in the Material
* is granted to  or  conferred  upon  you,  either   expressly,  by implication,
* inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
* property rights must be express and approved by Intel in writing.
*
* Unless otherwise agreed by Intel in writing,  you may not remove or alter this
* notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
* suppliers or licensors in any way.
*******************************************************************************/

/*
*  Content:
*  Intel MKL RCI (P)FGMRES ((Preconditioned) Flexible Generalized Minimal
*                                                       RESidual method) example
********************************************************************************/

/*---------------------------------------------------------------------------
*  Example program for solving non-symmetric indefinite system of equations
*  Fully advanced case: full functionality of RCI FGMRES solver is exploited
*---------------------------------------------------------------------------*/
#ifndef FGMRES_HPP
#define FGMRES_HPP 1
#include <stdio.h>
#include "mkl.h"
#include <cstdio>
#include <vector>
//#include"mkl.h"
//#include <mkl_blas.h>
//#include <mkl_lapacke.h>
//#include<mkl_types.h>
//#include "mkl_rci.h"
//#include "mkl_blas.h"
//#include "mkl_spblas.h"
//#include "mkl_service.h"
#include "fmm.hpp"
#include <Eigen/Dense>
#include "Eigen/Core"

using std::vector;
using namespace Eigen;
using namespace std;
namespace fmm_3d {
	int gmres(vector<zone>& Zone, GValue& gvalue, BEMGeometry& bemGeo, M2Ematrix* ME, E2Lmatrix* EL, math_pro& math_p, Multipole& op, double rhs[], int N, int numzone,int& interation)
	{

		const MKL_INT size = 300;

		/*---------------------------------------------------------------------------
		* Allocate storage for the ?par parameters and the solution/rhs/residual vectors
		*---------------------------------------------------------------------------*/
		MKL_INT ipar[size], * ipiv, info;
		double dpar[size];
		double* tmp;
		double* computed_solution, * b;
		/*---------------------------------------------------------------------------
		* Some additional variables to use with the RCI (P)FGMRES solver
		*---------------------------------------------------------------------------*/
		MKL_INT itercount;
		MKL_INT RCI_request, i, ivar, isw;
		double dvar;
		char cvar;
		double bnrm2;
		const int M = 250;
		printf("--------------------------------------------------------\n");
		printf("The FULLY ADVANCED example of usage of RCI FGMRES solver\n");
		printf("   to solve a non-symmetric indefinite non-degenerate\n");
		printf("          algebraic system of linear equations\n");
		printf("--------------------------------------------------------\n\n");
		/*---------------------------------------------------------------------------
		* Initialize variables and the right hand side through matrix-vector product
		*---------------------------------------------------------------------------*/
		ivar = N;
		cvar = 'T';
		i = 1;
		isw = 1;
		/*dcopy(&ivar, rhs, &i, b, &i);*/
		bnrm2 = dnrm2(&ivar, rhs, &i); cout << "oook" << endl;
		/*---------------------------------------------------------------------------
		* Initialize the initial guess
		*---------------------------------------------------------------------------*/
		//size_t a = N * (2 * M + 1) + (M * (M + 9)) / 2 + 1;
	 //   cout << "到了a" << a << endl;
		tmp = new double[N * (2 * M + 1) + (M * (M + 9)) / 2 + 1];
		cout << "开辟空间" << endl;
		computed_solution = new double[N];
		b = new double[N];
		ipiv = new MKL_INT[N];
		for (i = 0; i < N; i++)
		{
			computed_solution[i] = 0.0;
		}

		/*---------------------------------------------------------------------------
		* Initialize the solver
		*---------------------------------------------------------------------------*/
		dfgmres_init(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp);
		if (RCI_request != 0)
			goto FAILED;
		/*---------------------------------------------------------------------------
		* Set the desired parameters:
		* do the restart after 2 iterations
		* LOGICAL parameters:
		* do not do the stopping test for the maximal number of iterations
		* do the Preconditioned iterations of FGMRES method
		* DOUBLE PRECISION parameters
		* set the relative tolerance to 1.0D-3 instead of default value 1.0D-6
		*---------------------------------------------------------------------------*/
		//ipar[5] = 100000000000;//最大迭代次数
		//ipar[4] = ipar[5];
		ipar[14] = M;
		ipar[7] = 0;
		ipar[8] = 0;
		ipar[10] = 1;    //if ipar[10] is zero, the dfgmres routine runs the non - preconditioned version of the FGMRES method
		ipar[11] = 0;
		dpar[0] = gvalue.tol;
		/*---------------------------------------------------------------------------
		* Check the correctness and consistency of the newly set parameters
		*---------------------------------------------------------------------------*/
		//dfgmres_check(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp);

		if (RCI_request != 0)
			goto FAILED;
		/*---------------------------------------------------------------------------
		* Print the info about the RCI FGMRES method
		*---------------------------------------------------------------------------*/
		printf("Some info about the current run of RCI FGMRES method:\n\n");
		if (ipar[7])
		{
			printf("As ipar[7]=%d, the automatic test for the maximal number of ", ipar[7]);
			printf("iterations will be\nperformed\n");
		}
		else
		{
			printf("As ipar[7]=%d, the automatic test for the maximal number of ", ipar[7]);
			printf("iterations will be\nskipped\n");
		}
		printf("+++\n");
		if (ipar[8])
		{
			printf("As ipar[8]=%d, the automatic residual test will be performed\n", ipar[8]);
		}
		else
		{
			printf("As ipar[8]=%d, the automatic residual test will be skipped\n", ipar[8]);
		}
		printf("+++\n");
		if (ipar[9])
		{
			printf("As ipar[9]=%d, the user-defined stopping test will be ", ipar[9]);
			printf("requested via\nRCI_request=2\n");
		}
		else
		{
			printf("As ipar[9]=%d, the user-defined stopping test will not be ", ipar[9]);
			printf("requested, thus,\nRCI_request will not take the value 2\n");
		}
		printf("+++\n");
		if (ipar[10])
		{
			printf("As ipar[10]=%d, the Preconditioned FGMRES iterations will be ", ipar[10]);
			printf("performed, thus,\nthe preconditioner action will be requested via ");
			printf("RCI_request=3\n");
		}
		else
		{
			printf("As ipar[10]=%d, the Preconditioned FGMRES iterations will not ", ipar[10]);
			printf("be performed,\nthus, RCI_request will not take the value 3\n");
		}
		printf("+++\n");
		if (ipar[11])
		{
			printf("As ipar[11]=%d, the automatic test for the norm of the next ", ipar[11]);
			printf("generated vector is\nnot equal to zero up to rounding and ");
			printf("computational errors will be performed,\nthus, RCI_request will not ");
			printf("take the value 4\n");
		}
		else
		{
			printf("As ipar[11]=%d, the automatic test for the norm of the next ", ipar[11]);
			printf("generated vector is\nnot equal to zero up to rounding and ");
			printf("computational errors will be skipped,\nthus, the user-defined test ");
			printf("will be requested via RCI_request=4\n");
		}
		printf("+++\n\n");
		/*---------------------------------------------------------------------------
		* Compute the solution by RCI (P)FGMRES solver with preconditioning
		* Reverse Communication starts here
		*---------------------------------------------------------------------------*/

	ONE:dfgmres(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp);

		/*---------------------------------------------------------------------------
		* If RCI_request=0, then the solution was found with the required precision
		*---------------------------------------------------------------------------*/
		if (RCI_request == 0)
			goto COMPLETE;
		/*---------------------------------------------------------------------------
		* If RCI_request=1, then compute the vector A*tmp[ipar[21]-1]
		* and put the result in vector tmp[ipar[22]-1]
		*---------------------------------------------------------------------------
		* NOTE that ipar[21] and ipar[22] contain FORTRAN style addresses,
		* therefore, in C code it is required to subtract 1 from them to get C style
		* addresses
		*---------------------------------------------------------------------------*/
		if (RCI_request == 1)
		{
		/*	for (int i = 0; i < N; ++i)
			{
				cout << "赋值赋错了？" << bemGeo.Panels[i].q << endl;
			}*/

			int count = 0;

			for (int i = 0; i < N; ++i)
			{
	

				if (bemGeo.Panels[gvalue.iElem[i]].bctype == 1)
				{
					bemGeo.Panels[gvalue.iElem[i]].phi = tmp[ipar[21] - 1 + count];
					count++;
					
				}
				else if (bemGeo.Panels[gvalue.iElem[i]].bctype == 2)
				{
					bemGeo.Panels[gvalue.iElem[i]].q = tmp[ipar[21] - 1 + count];
					count++;

				}
				else if (bemGeo.Panels[gvalue.iElem[i]].bctype == 3)
				{
					bemGeo.Panels[gvalue.iElem[i]].q = tmp[ipar[21] - 1 + count];
					count++;
					bemGeo.Panels[gvalue.iElem[i]].phi = tmp[ipar[21] - 1 + count];
				
					count++;

				}
				else if (bemGeo.Panels[gvalue.iElem[i]].bctype == 4)
				{
					bemGeo.Panels[gvalue.iElem[i]].phi = bemGeo.Panels[gvalue.iElem[i]].dualPanel->phi;
					bemGeo.Panels[gvalue.iElem[i]].q = bemGeo.Panels[gvalue.iElem[i]].dualPanel->q;
					//cout << "单元编号" << gvalue.iElem[i] << endl;//存的编号重了
				}

				else
				{
					cout << "this is error!" << endl;
				}


			}

			//for (int i = 0; i < N; ++i)
			//{
			//	cout << "赋值赋错了？" << bemGeo.Panels[i].q <<"tmp"<< tmp[ipar[21] - 1 +i]<< endl;
			//}

			int b_loc = 0;
			for (int i = 0; i < N; ++i)
			{
				tmp[ipar[22] - 1 + i] = 0.0;

			}
			//这出问题了
			for (int j = 0; j < numzone; ++j)
			{
				int k = gvalue.turn[j] - 1;
				op.upward_pass( gvalue, bemGeo, Zone[k], ME, math_p, isw);
				op.downward_pass( gvalue, bemGeo, Zone[k], EL, math_p, &tmp[ipar[22] - 1], isw);
				//cout << "求解完的域" << k << endl;
			}


			//Eigen::MatrixXd rot2(N, N);
			//for (int i = 0; i < N; ++i)
			//{
			//	for (int j = 0; j < N; ++j)
			//	{
			//		rot2(i, j) = 0.0;
			//	}
			//}
			//



			if (isw == 1) {
				cout << "迭代次数" << isw << endl;

				
				for (int k = 0; k < numzone; ++k)
				{
					int j = gvalue.turn[k] - 1;
					for (int i = 0; i < Zone[j].id; ++i) {
						if (Zone[j].cell[i].isLeaf) {
							if (Zone[j].cell[i].pretype == 1)
							{
				
								info = LAPACKE_dgetrf(LAPACK_COL_MAJOR, Zone[j].cell[i].nElem, Zone[j].cell[i].nElem, &(gvalue.amat[Zone[j].cell[i].preLoc]), Zone[j].cell[i].nElem, &ipiv[Zone[j].cell[i].xLoc]);
								//for (int m = 0; m < Zone[j].cell[i].nElem; ++m)
								//{
								//	for (int n = 0; n < Zone[j].cell[i].nElem; ++n)
								//	{
								//		//cout << "行" << Zone[j].cell[i].xLoc + m << "列" << Zone[j].cell[i].xLoc + n << endl;
								//		rot2(Zone[j].cell[i].xLoc + m, Zone[j].cell[i].xLoc + n) = gvalue.amat[Zone[j].cell[i].preLoc + m * Zone[j].cell[i].nElem + n];
								//	}
								//}
							}
							if (Zone[j].cell[i].pretype == 2)
							{
					
								info = LAPACKE_dgetrf(LAPACK_COL_MAJOR, Zone[j].cell[i].newLoc, Zone[j].cell[i].newLoc, &(gvalue.amat[Zone[j].cell[i].preLoc]), Zone[j].cell[i].newLoc, &ipiv[Zone[j].cell[i].xLoc]);
								/*for (int m = 0; m < Zone[j].cell[i].newLoc; ++m)
								{
									for (int n = 0; n < Zone[j].cell[i].newLoc; ++n)
									{
										rot2(Zone[j].cell[i].xLoc + m, Zone[j].cell[i].xLoc + n) = gvalue.amat[Zone[j].cell[i].preLoc + m * Zone[j].cell[i].newLoc + n];
									}
								}*/
							}

						
						}
					
					}
				
				}
			

			}

			/*std::ofstream fout;
			fout.open("对角矩阵A.dat");
			for (int i = 0; i < N; ++i)
			{
				for (int j = 0; j < N; ++j)

				{
					fout << rot2(i,j) << endl;
				}


			}
			fout.close();

			fout.open("一维数组.dat");
			for (int i = 0; i < gvalue.preSize; ++i)
			{
			
					fout << gvalue.amat[i] << endl;
				


			}
			fout.close();*/

			isw = 2;
			interation++;
			goto ONE;
		}
		/*---------------------------------------------------------------------------
		* If RCI_request=2, then do the user-defined stopping test
		* The residual stopping test for the computed solution is performed here
		*---------------------------------------------------------------------------
		* NOTE: from this point vector b[N] is no longer containing the right-hand
		* side of the problem! It contains the current FGMRES approximation to the
		* solution. If you need to keep the right-hand side, save it in some other
		* vector before the call to dfgmres routine. Here we saved it in vector
		* rhs[N]. The vector b is used instead of rhs to preserve the
		* original right-hand side of the problem and guarantee the proper
		* restart of FGMRES method. Vector b will be altered when computing the
		* residual stopping criterion!
		*---------------------------------------------------------------------------*/
		if (RCI_request == 2)
		{
			/* Request to the dfgmres_get routine to put the solution into b[N] via ipar[12]
			   --------------------------------------------------------------------------------
			   WARNING: beware that the call to dfgmres_get routine with ipar[12]=0 at this
			   stage may destroy the convergence of the FGMRES method, therefore, only
			   advanced users should exploit this option with care */
			ipar[12] = 1;
			/* Get the current FGMRES solution in the vector b[N] */

			dfgmres_get(&ivar, computed_solution, b, &RCI_request, ipar, dpar, tmp, &itercount);


			/* Compute the current true residual via MKL (Sparse) BLAS routines */
			//mkl_dcsrgemv (&cvar, &ivar, A, ia, ja, b, residual);
			//dvar = -1.0E0;
			//i = 1;
			//daxpy (&ivar, &dvar, rhs, &i, residual, &i);
			//dvar = dnrm2 (&ivar, residual, &i);
			printf("iteration no. = %d, relative residual = %lf\n", itercount, dpar[4] / bnrm2);
			if (dpar[4] / bnrm2 < dpar[0]) {
				goto COMPLETE;
			}
			else {
				goto ONE;
			}
		}
		/*---------------------------------------------------------------------------
		* If RCI_request=3, then apply the preconditioner on the vector
		* tmp[ipar[21]-1] and put the result in vector tmp[ipar[22]-1]
		*---------------------------------------------------------------------------
		* NOTE that ipar[21] and ipar[22] contain FORTRAN style addresses,
		* therefore, in C code it is required to subtract 1 from them to get C style
		* addresses
		*---------------------------------------------------------------------------*/
		if (RCI_request == 3)
		{
			//add your code here ....
			for (int i = 0; i < N; ++i) {
				tmp[ipar[22] - 1 + i] = tmp[ipar[21] - 1 + i];
			}
//开始预处理
			int preLoc = 0;
			for (int k = 0; k < numzone; ++k)
			{
				int j = gvalue.turn[k] - 1;
				for (int i = 0; i < Zone[j].id; ++i) {
					if (Zone[j].cell[i].isLeaf) {
						if (Zone[j].cell[i].pretype == 1)
						{
						LAPACKE_dgetrs(LAPACK_COL_MAJOR, cvar, Zone[j].cell[i].nElem, 1, &(gvalue.amat[Zone[j].cell[i].preLoc]),
								Zone[j].cell[i].nElem, &ipiv[Zone[j].cell[i].xLoc], &tmp[ipar[22] - 1 + Zone[j].cell[i].xLoc], Zone[j].cell[i].nElem);
						//cout << "1111111我是这里出了问题" << endl;
						}
						if (Zone[j].cell[i].pretype == 2)//这个也不对，因为是2，3，2，3
						{
							LAPACKE_dgetrs(LAPACK_COL_MAJOR, cvar, Zone[j].cell[i].newLoc, 1, &(gvalue.amat[Zone[j].cell[i].preLoc]),
								Zone[j].cell[i].newLoc, &ipiv[Zone[j].cell[i].xLoc], &tmp[ipar[22] - 1 + Zone[j].cell[i].xLoc], Zone[j].cell[i].newLoc);
							//cout << "我是这里出了问题" << endl;
						}
					}
				}
			
			}


			goto ONE;
		}
		/*---------------------------------------------------------------------------
		* If RCI_request=4, then check if the norm of the next generated vector is
		* not zero up to rounding and computational errors. The norm is contained
		* in dpar[6] parameter
		*---------------------------------------------------------------------------*/
		if (RCI_request == 4)
		{
			if (dpar[6] < 1.0E-12) {
				goto COMPLETE;
			}
			else {
				goto ONE;
			}
		}
		/*---------------------------------------------------------------------------
		* If RCI_request=anything else, then dfgmres subroutine failed
		* to compute the solution vector: computed_solution[N]
		*---------------------------------------------------------------------------*/
		else
		{
			goto FAILED;
		}
		/*---------------------------------------------------------------------------
		* Reverse Communication ends here
		* Get the current iteration number and the FGMRES solution (DO NOT FORGET to
		* call dfgmres_get routine as computed_solution is still containing
		* the initial guess!). Request to dfgmres_get to put the solution
		* into vector computed_solution[N] via ipar[12]
		*---------------------------------------------------------------------------*/
	COMPLETE:ipar[12] = 1;
		dfgmres_get(&ivar, computed_solution, rhs, &RCI_request, ipar, dpar, tmp, &itercount);

		printf(" The system has been solved \n");
		delete[] computed_solution, tmp, b;
		MKL_Free_Buffers();
		return 0;

		/*-------------------------------------------------------------------------*/
		/* Release internal MKL memory that might be used for computations         */
		/* NOTE: It is important to call the routine below to avoid memory leaks   */
		/* unless you disable MKL Memory Manager                                   */
		/*-------------------------------------------------------------------------*/
	FAILED:printf("\nThis example FAILED as the solver has returned the ERROR code %d", RCI_request);
		delete[] computed_solution, tmp, b;
		MKL_Free_Buffers();
		return 1;
	}
}
#endif