C++ ComplexMatrix类代码示例

void NRLib::ComputeEigenVectorsComplex(ComplexMatrix        & A,
                                       ComplexVector        & eigen_values,
                                       ComplexMatrix        & eigen_vectors)
{
  ComplexMatrix dummy_mat(A.numRows(), A.numCols());
  flens::ev(false, true, A, eigen_values, dummy_mat, eigen_vectors);
}

void NRLib::Adjoint(const ComplexMatrix & in,
                    ComplexMatrix       & out)
{
  int m = out.numRows();
  int n = out.numCols();

  for (int i=0 ; i < m ; i++) {
    for (int j=0 ; j < n ; j++) {
      out(i,j) = std::conj(in(j,i));
    }
  }
}

ComplexMatrix MatrixSqrtHermitean (const ComplexMatrix& A)
//
// Returns the matrix square root of a hermitean matrix A.
// The argument matrix is not checked for hermiticity !
//
//
{  
    // get the dimension information
    int lo = A.Clo(),
        hi = A.Chi();

    // columns and rows must have the same range
    if (A.Rlo() != lo || A.Rhi() != hi) 
        Matpack.Error(Mat::UnspecifiedError,"MatrixSqrtHermitean: hermitean matrix must be square"); 
    
    Matrix z(lo,hi,lo,hi), 
           zr(lo,hi,lo,hi),
           zi(lo,hi,lo,hi);

    Vector d(lo,hi);

    int i,j;

    // create packed form of hermitean matrix A in matrix z
    for (i = lo; i <= hi; i++) {        
        for (j = lo; j < i; j++) {
            z[j][i]=imag(A[i][j]);
            z[i][j]=real(A[i][j]);
        }
        z[i][i]=real(A[i][i]);
    }
    
    // diagonalize z
    EigenSystemHermitean(z,d,zr,zi,false,30);

    // combine real and imaginary parts of eigenvectors (zr,zi) in complex matrix
    ComplexMatrix U(zr,zi);
  
    ComplexVector dc(lo,hi);

    for (i = lo; i <= hi; i++) {        
        if (d[i] < 0) 
          dc[i]=complex<double>(0,sqrt(-d[i]));
        else 
          dc[i]=complex<double>(sqrt(d[i]),0);
    }

    // calculate matrix square root
    U = U * Diagonal(dc) * U.Hermitean();

    return U.Value();    
}

void NRLib::WriteComplexMatrix(const std::string   & header,
                               const ComplexMatrix & C)
{
  int m = C.numRows();
  int n = C.numCols();
  LogKit::LogFormatted(LogKit::Error,"\n"+header+"\n");
  for (int i=0; i < m ; i++) {
    for (int j=0; j < n ; j++) {
      LogKit::LogFormatted(LogKit::Error,"(%12.8f, %12.8f) ",C(i,j).real(),C(i,j).imag());
    }
    LogKit::LogFormatted(LogKit::Error,"\n");
  }
  LogKit::LogFormatted(LogKit::Error,"\n");
}

static void transform(ComplexMatrix &inout, const std::vector<int> &loop_dims_select, unsigned fftw_flags, bool forward)
{
	if (os::disable_SSE_for_FFTW()) {
		fftw_flags |= FFTW_UNALIGNED; // see os.h
	}

	int num_dims, num_loop_dims;
	fftw_iodim dims[16], loop_dims[16];
	make_iodims(inout.getShape(), loop_dims_select, num_dims, dims, num_loop_dims, loop_dims);

	ComplexMatrix::accessor inout_acc(inout);
	double *re = inout_acc.ptr_real();
	double *im = inout_acc.ptr_imag();

	if (!forward) {
		std::swap(re, im);
	}

	fftw_plan plan = fftw_plan_guru_split_dft(
		num_dims, dims, 
		num_loop_dims, loop_dims,
		re, im, // in
		re, im, // out
		fftw_flags
	);
	assert(plan);

	fftw_execute_split_dft(plan, re, im, re, im);

	fftw_destroy_plan(plan);
}

DoubleMatrix mult(DoubleMatrix& m2, ComplexMatrix& m1)
{
    //  Check dimensions
    unsigned int m1_nRows = m1.numRows();
    unsigned int m2_nRows = m2.numRows();
    unsigned int m1_nColumns = m1.numCols();
    unsigned int m2_nColumns = m2.numCols();

    if (m1.size() == 0)
    {
        return real(m1);
    }

    if (m2.size() == 0)
    {
        return m2;
    }

    DoubleMatrix result(m1_nRows, m2_nColumns);
    if (m1_nColumns == m2_nRows)
    {
        for (unsigned int row = 0; row < result.numRows(); row++)
        {
            for (unsigned int col = 0; col < m2_nColumns; col++)
            {
                double sum = 0.0;
                for (unsigned int k = 0; k < m1_nColumns; k++)
                {
                    sum = sum + (real(m1[row][k]) * m2[k][col]);
                }
                result[row][col] = sum;
            }
        }
        return result;
    }

    if (m1_nRows == m2_nColumns)
    {
        return mult(m2, m1);
    }

    throw ("Incompatible matrix operands to multiply");
}

ComplexMatrix subtract(ComplexMatrix& x, ComplexMatrix& y)
{
    if(sameDimensions(x,y))
    {
        ComplexMatrix result(x.RSize(), x.CSize());

        for (int i = 0; i < x.RSize(); i++)
        {
            for (int j = 0; j < x.CSize(); j++)
            {
                result(i, j) = x(i, j) - y(i, j);
            }
        }
        return result;
    }
    else
    {
       throw ("Matrices must be the same dimension to perform subtraction");
    }
}

// Matrix multiplication definition
ComplexMatrix operator*(const ComplexMatrix &A, const ComplexMatrix& B)
{
    assert(A.IsSquare());

    int numRows = A.GetNumberOfRows();
    int numCols = A.GetNumberOfCols();
    ComplexMatrix C(numRows, numCols);

    for (int i  = 0; i < numRows; i++)
    {
        for (int j = 0; j < numRows; j++)
        {
            for (int k = 0; k < numRows; k++)
            {
                C.mMemory[i][j] = C.mMemory[i][j] + A.mMemory[i][k] * B.mMemory[k][j];
            }
        }
    }

    return C;
}

void run_test(const Matrix &mat, int k, int m)
{
    DenseGenMatProd<double> op(mat);
    GenEigsSolver<double, SelectionRule, DenseGenMatProd<double>> eigs(&op, k, m);
    eigs.init();
    int nconv = eigs.compute();
    int niter = eigs.num_iterations();
    int nops = eigs.num_operations();

    REQUIRE( nconv > 0 );

    ComplexVector evals = eigs.eigenvalues();
    ComplexMatrix evecs = eigs.eigenvectors();

    ComplexMatrix err = mat * evecs - evecs * evals.asDiagonal();

    INFO( "nconv = " << nconv );
    INFO( "niter = " << niter );
    INFO( "nops = " << nops );
    INFO( "||AU - UD||_inf = " << err.array().abs().maxCoeff() );
    REQUIRE( err.array().abs().maxCoeff() == Approx(0.0) );
}

ls::ComplexMatrix getComplexMatrixFromString(const string& textMatrix)
{
	ComplexMatrix mat;

    //Parse the matrix
    vector<string> rows = splitString(textMatrix, "\n");
    for(int row = 0; row < rows.size(); row++)
    {
        vector<string> values = splitString(rows[row], " \t");
        for(int col = 0; col < values.size(); col++)
        {
        	if(!mat.size())
            {
                mat.resize(rows.size(), values.size());
            }

            std::complex<double> val = toComplex(values[col]);
            mat(row, col).Real = real(val);
			mat(row, col).Imag = imag(val);
        }
    }
	return mat;
}

void run_test(const MatType& mat, int k, int m, double sigma, bool allow_fail = false)
{
    typename OpTypeTrait<MatType>::OpType op(mat);
    GenEigsRealShiftSolver<double, SelectionRule, typename OpTypeTrait<MatType>::OpType>
        eigs(&op, k, m, sigma);
    eigs.init();
    int nconv = eigs.compute();
    int niter = eigs.num_iterations();
    int nops  = eigs.num_operations();

    if(allow_fail)
    {
        if( eigs.info() != SUCCESSFUL )
        {
            WARN( "FAILED on this test" );
            std::cout << "nconv = " << nconv << std::endl;
            std::cout << "niter = " << niter << std::endl;
            std::cout << "nops  = " << nops  << std::endl;
            return;
        }
    } else {
        INFO( "nconv = " << nconv );
        INFO( "niter = " << niter );
        INFO( "nops  = " << nops );
        REQUIRE( eigs.info() == SUCCESSFUL );
    }

    ComplexVector evals = eigs.eigenvalues();
    ComplexMatrix evecs = eigs.eigenvectors();

    ComplexMatrix resid = mat * evecs - evecs * evals.asDiagonal();
    const double err = resid.array().abs().maxCoeff();

    INFO( "||AU - UD||_inf = " << err );
    REQUIRE( err == Approx(0.0) );
}

//---------------------------------------------------------
void NRLib::CholeskySolveComplex(ComplexMatrix & A,
                                 ComplexMatrix & B) // This is also where the solution is stored
//---------------------------------------------------------
{
  NRLib::ComplexMatrix::TriangularView::HermitianView H = A.lower().hermitian();

  int info = flens::posv(H, B);

  if (info != 0) {
    std::ostringstream oss;
    if (info < 0) {
      oss << "Internal FLENS/Lapack error: Error in argument " << -info
          << " of posv call.";
    }
    else {  // info > 0
      oss << "Error in Cholesky: The leading minor of order " << info
          << " is not positive definite.";
    }
    throw Exception(oss.str());
  }
}

TEST(AlphaBetaTest, FourSitesNoDisorder) {
	Parameters pars;
	pars.nmax = 3;
	pars.e0 = pars.t0 = pars.d0 = 1;
	pars.e0MaxDisorder = pars.t0MaxDisorder = pars.d0MaxDisorder = 0.0;
	pars.e0seed = pars.t0seed = pars.d0seed = 1;

	AlphaBeta ab(pars);
	int nsum = 1;
	complex_mkl z = {1.0, 0.1};
	ComplexMatrix alpha;
	ComplexMatrix beta;
	ab.FillAlphaBetaMatrix(nsum, z, alpha, beta);
	// testing the dimensions
	EXPECT_EQ(alpha.GetRows(),0);
	EXPECT_EQ(alpha.GetCols(),0);
	EXPECT_EQ(beta.GetRows(),1);
	EXPECT_EQ(beta.GetCols(),1);
	//testing matrix element
	dcomplex x = pars.t0/(convertToDcomplex(z)-pars.e0-pars.e0-pars.d0);
	EXPECT_DOUBLE_EQ(beta(0,0).real,x.real());
	EXPECT_DOUBLE_EQ(beta(0,0).imag,x.imag());

	nsum=2;
	ab.FillAlphaBetaMatrix(nsum, z, alpha, beta);
	// testing the dimensions
	EXPECT_EQ(alpha.GetRows(),1);
	EXPECT_EQ(alpha.GetCols(),1);
	EXPECT_EQ(beta.GetRows(),1);
	EXPECT_EQ(beta.GetCols(),2);
	//testing matrix element
	x = pars.t0/(convertToDcomplex(z)-pars.e0-pars.e0);
	EXPECT_DOUBLE_EQ(alpha(0,0).real,x.real());
	EXPECT_DOUBLE_EQ(alpha(0,0).imag,x.imag());
	EXPECT_DOUBLE_EQ(beta(0,0).real,x.real());
	EXPECT_DOUBLE_EQ(beta(0,0).imag,x.imag());
	EXPECT_DOUBLE_EQ(beta(0,1).real,x.real());
	EXPECT_DOUBLE_EQ(beta(0,1).imag,x.imag());

	nsum =3;
	ab.FillAlphaBetaMatrix(nsum, z, alpha, beta);
	// testing the dimensions
	EXPECT_EQ(alpha.GetRows(),2);
	EXPECT_EQ(alpha.GetCols(),1);
	EXPECT_EQ(beta.GetRows(),2);
	EXPECT_EQ(beta.GetCols(),1);
	//testing matrix element
	x = pars.t0/(convertToDcomplex(z)-pars.e0-pars.e0);
	EXPECT_DOUBLE_EQ(alpha(0,0).real,x.real());
	EXPECT_DOUBLE_EQ(alpha(0,0).imag,x.imag());
	x = pars.t0/(convertToDcomplex(z)-pars.e0-pars.e0-pars.d0);
	EXPECT_DOUBLE_EQ(alpha(1,0).real,x.real());
	EXPECT_DOUBLE_EQ(alpha(1,0).imag,x.imag());
	x = pars.t0/(convertToDcomplex(z)-pars.e0-pars.e0);
	EXPECT_DOUBLE_EQ(beta(0,0).real,x.real());
	EXPECT_DOUBLE_EQ(beta(0,0).imag,x.imag());
	x = pars.t0/(convertToDcomplex(z)-pars.e0-pars.e0-pars.d0);
	EXPECT_DOUBLE_EQ(beta(1,0).real,x.real());
	EXPECT_DOUBLE_EQ(beta(1,0).imag,x.imag());

	nsum =4;
	ab.FillAlphaBetaMatrix(nsum, z, alpha, beta);
	// testing the dimensions
	EXPECT_EQ(alpha.GetRows(),1);
	EXPECT_EQ(alpha.GetCols(),2);
	EXPECT_EQ(beta.GetRows(),1);
	EXPECT_EQ(beta.GetCols(),1);
	//testing matrix element
	x = pars.t0/(convertToDcomplex(z)-pars.e0-pars.e0);
	EXPECT_DOUBLE_EQ(alpha(0,0).real,x.real());
	EXPECT_DOUBLE_EQ(alpha(0,0).imag,x.imag());
	EXPECT_DOUBLE_EQ(alpha(0,1).real,x.real());
	EXPECT_DOUBLE_EQ(alpha(0,1).imag,x.imag());
	EXPECT_DOUBLE_EQ(beta(0,0).real,x.real());
	EXPECT_DOUBLE_EQ(beta(0,0).imag,x.imag());


	nsum = 5;
	ab.FillAlphaBetaMatrix(nsum, z, alpha, beta);
	// testing the dimensions
	EXPECT_EQ(alpha.GetRows(),1);
	EXPECT_EQ(alpha.GetCols(),1);
	EXPECT_EQ(beta.GetRows(),0);
	EXPECT_EQ(beta.GetCols(),0);
	//testing matrix element
	x = pars.t0/(convertToDcomplex(z)-pars.e0-pars.e0-pars.d0);
	EXPECT_DOUBLE_EQ(alpha(0,0).real,x.real());
	EXPECT_DOUBLE_EQ(alpha(0,0).imag,x.imag());


}

//From Fransk CSharp code... a bug was fixed in this code..
bool sameDimensions(ComplexMatrix& x, ComplexMatrix& y)
{
    return ((x.RSize() == y.RSize()) && (x.CSize() == y.CSize())) ? true : false;
}

int main (void)
{
  MpImage image(size,size);		     // allocate size x size raster image

  MpImage backimg;			                 // read background image
  if (back_image) {
    if ( ! backimg.ReadAnyFile(back_image) )
      Matpack.Error("Can't read \"%s\"",back_image);
  }

  // read potential matrix
  ifstream pot_stream("pot");
  if ( ! pot_stream) Matpack.Error("Can't open pot data file \"%s\"","pot");
  pot_stream >> pot;
  pot_stream.close();

  for (int i = 0; i <= 175; i += 5) {  			   // loop through frames

    fprintf(stdout,"\rframe %d...",i); fflush(stdout);         // what's going on

    // read complex wave function
    char file_name[200], pipe_name[200], output_name[200];
    sprintf(file_name,"psi%d.gz",i);		           // generate input data file name
    sprintf(pipe_name,"gunzip -c %s > psi.dat",file_name); // decompress data file 
    sprintf(output_name,"psi%03d.jpg",i);                  // generate output file name

    system(pipe_name);
    ifstream data_stream("psi.dat");
    if ( ! data_stream) Matpack.Error("Can't open data file psi.dat");
    data_stream >> psi;
    data_stream.close();
    unlink("psi.dat");

    // consistency checks
    if ( pot.Rlo() != psi.Rlo() || pot.Rhi() != psi.Rhi() || 
	 pot.Clo() != psi.Clo() || pot.Chi() != psi.Chi() )
      Matpack.Error("non-conformant index range of potential and wave function");

    Scene scene;				   // define scene for 3d drawing

    // create surface
    double u0 = psi.Clo(), u1 = psi.Chi(), 
           v0 = psi.Rlo(), v1 = psi.Rhi();
    int    nu = psi.Chi()-psi.Clo(), nv = psi.Rhi()-psi.Rlo(),  
           su = 0, sv = 0,  periodic = 0;
    ParametricSurface(scene, height_fcn, u0,u1,nu,su, v0,v1,nv,sv, 
		      periodic, Identity, color_fcn);

    scene.BoxRatios(1,z_aspect,1);               // scale scene to box with x:y:z

    float dist = 1.6;			                    // distance parameter
    scene.Look(Vector3D( 1.1*dist, 1.5*dist, 1.4*dist), // camera position vector
	       Vector3D(-1.1*dist,-1.65*dist,-1.4*dist), // look direction vector
	       FieldOfView(45), 		 // field of view angle in degree
	       0);		             // "twist your head" angle in degree

    scene.SetGlobalShading(Facet::Gouraud);   
                                       // shading (None, Flat, Gouraud, Phong...)
    scene.SetHiddenSurface(preview ? Scene::DepthSort : Scene::TopoSort);

    // edgeline drawing option (None,...)
    scene.SetGlobalEdgeLines(show_grid ? Facet::Individual : None);

    // Setting per-vertex coloring is the best choice to get smoothly changing
    // colors on the surface. This requires Gouraud or Phong shading to be
    // switched on. If you use per-facet coloring then the facet edges are still
    // slightly visible ("Scene::PerFacet")
    scene.SetColoring(Scene::PerVertex);

    // Define ambient light, RGB color, and specular light and exponent
    // If per-vertex coloring or per-facet coloring is set, then only the ambient
    // light factor and the specular light factor and exponent are in effect, and
    // the diffuse color is ignored (the vertex or facet color is used instead)
    Material material(0.6, ColorF(1.,0.8,0.), 0.3,1.8); 
    scene.SetFacetMaterial(material);		    // set all facets to material

    scene.Open(image);			        // direct drawing to raster image
    if (back_image)
      image.InsertTiled(backimg);
    else
      scene.SetBackground(back_color); 	        // draw background with RGB color
    scene.SetColor(edge_color);		           // define color for edge lines
    scene.Show();			  		    // render the surface
    scene.Close();				       // call always after close

    image.WriteJpegFile(output_name,jpeg_quality,jpeg_smooth);// write JPEG image

    // write GIF image, if you prefer that
    //image.WriteGifFile(output_name);     // write image as GIF

  }

  cout << "\nfinished" << endl;
}