psi::extrema::math_tools Class Reference

A general mathematics routine class. More...

#include <math_tools.h>

Inheritance diagram for psi::extrema::math_tools:

Protected Member Functions
	math_tools ()
	Default constructor, does nothing.
	~math_tools ()
	Default destructor, does nothing.
double *	newton_step (int dim, double **_Hi, double *g)
	Computes the Newton-Raphson optimization step.
double **	update_bfgs (int dim, double *_var_dif, double *grad_dif, double **_Hi_old)
	Performs bfgs update on inverse of hessian matrix.
double **	update_ms (int dim, double *_var_dif, double *_grad_dif, double **_Hi_old)
	Performs update on inverse of hessian matrix attributed to Murtah and Sargent among many others.
double **	rep_reduce (char *label, double **rep_matrix, int reps)
	Reduces a set of representation vectors.
double **	rep_project (char *label, int dim_vec, double **result_vecs, int *irrep_proj)
	Projects out irrep components from a reducible representation.
double **	orthogonalize (int nvecs, int dimvecs, double **vecs, int normalize, double norm_tol, int *nindep)
	Othogonalizes a set of vectors.

---

Detailed Description

A general mathematics routine class.

General implementations of basic mathematical algorithms. The only dependency other than standard libraries is the PSI 3.0 library libciomr.

Definition at line 22 of file math_tools.h.

---

Member Function Documentation

double * math_tools::newton_step	(	int	dim,
		double **	_Hi,
		double *	g
	)			[protected]

Computes the Newton-Raphson optimization step.

Parameters:

	dim	number of variables
	**_Hi	inverse hessian matrix
	*g	gradient vector

Returns:

*_s displacement vector

Definition at line 25 of file math_tools.cc.

References init_array().

Referenced by psi::extrema::zmat::newton_step().

                                                                  {
    
    int i, j;
    double *_s;
        
    _s = init_array(dim);
    
    for(i=0;i<dim;++i) 
        for(j=0;j<dim;++j) 
            _s[i] += -1.0 * _Hi[i][j] * g[j];

    return _s;
}

double ** math_tools::update_bfgs	(	int	dim,
		double *	_var_dif,
		double *	_grad_dif,
		double **	_Hi_old
	)			[protected]

Performs bfgs update on inverse of hessian matrix.

Parameters:

	dim	dimension of hessian inverse
	*_var_dif	difference of current and previous variable
	*grad_dif	difference of current and previous gradient
	**_Hi_old	inverse hessian from previous iteration

Returns:

**_Hi_new updated hessian inverse

Definition at line 52 of file math_tools.cc.

References free_matrix(), init_matrix(), and mmult().

Referenced by psi::extrema::coord_base::update_Hi().

                                                                        {

    int i,j;
    double **temp_mat, **_Hi_new, **mat1, **mat2, **mat3, num1=0.0, num2=0.0;
    _Hi_new = init_matrix(dim,dim);
    temp_mat = init_matrix(1,dim);
    mat1 = init_matrix(dim,dim);
    mat2 = init_matrix(dim,dim);
    mat3 = init_matrix(dim,dim);


    mmult(&_grad_dif,0,_Hi_old,0,temp_mat,0,1,dim,dim,0);
    for(i=0;i<dim;++i) 
        num1 += temp_mat[0][i]*_grad_dif[i];
 
    for(i=0;i<dim;++i)
        num2 += _var_dif[i]*_grad_dif[i];

    for(i=0;i<dim;++i)
        for(j=0;j<dim;++j)
            mat1[i][j] = _var_dif[i]*_var_dif[j];

    mmult(&_grad_dif,0,_Hi_old,0,temp_mat,0,1,dim,dim,0);
    for(i=0;i<dim;++i)
        for(j=0;j<dim;++j)
            mat2[i][j] = _var_dif[i]*temp_mat[0][j];

    free_matrix(temp_mat,1);
    temp_mat = init_matrix(dim,dim);
    for(i=0;i<dim;++i) 
        for(j=0;j<dim;++j)
            temp_mat[i][j] = _grad_dif[i]*_var_dif[j];
    mmult(_Hi_old,0,temp_mat,0,mat3,0,dim,dim,dim,0);

    for(i=0;i<dim;++i)
        for(j=0;j<dim;++j)
            _Hi_new[i][j] = _Hi_old[i][j] + (1 + num1/num2)*(mat1[i][j]/num2)
                - mat2[i][j]/num2 - mat3[i][j]/num2;
    
    free_matrix(temp_mat,dim);
    free_matrix(mat1,dim);
    free_matrix(mat2,dim);
    free_matrix(mat3,dim);

    return _Hi_new;
}

double ** math_tools::update_ms	(	int	dim,
		double *	_var_dif,
		double *	_grad_dif,
		double **	_Hi_old
	)			[protected]

Performs update on inverse of hessian matrix attributed to Murtah and Sargent among many others.

Parameters:

	dim	dimension of hessian inverse
	*_var_dif	difference of current and previous variable
	*_grad_dif	difference of current and previous gradient
	**_Hi_old	inverse hessian from previous iteration

Returns:

**_Hi_new = updated hessian inverse

Definition at line 114 of file math_tools.cc.

References free_matrix(), init_array(), and init_matrix().

Referenced by psi::extrema::coord_base::update_Hi().

                                                                        {

    int i, j;
    double div, *temp_arr, **temp_mat, **_Hi_new;

    /*allocate memory*/
    temp_arr = init_array(dim);
    temp_mat = init_matrix(dim,dim);
    _Hi_new = init_matrix(dim,dim);
    
    for(i=0;i<dim;++i) {
        for(j=0;j<dim;++j) {
            temp_arr[i] += _Hi_old[i][j] * _grad_dif[j];
        }
    }
    
    for(i=0;i<dim;++i) {
        temp_arr[i] = 0.0;
        for(j=0;j<dim;++j)
            temp_arr[i] += _Hi_old[i][j] * _grad_dif[j]; 
    }
    for(i=0;i<dim;++i)
        temp_arr[i] = _var_dif[i] - temp_arr[i];
    
    for(i=0;i<dim;++i)
        for(j=0;j<dim;++j)
            temp_mat[i][j] = temp_arr[i]*temp_arr[j];
    
    div = 0.0;
    for(i=0;i<dim;++i)
        div += temp_arr[i] * _grad_dif[i];
    
    for(i=0;i<dim;++i)
        for(j=0;j<dim;++j)
            _Hi_new[i][j] = _Hi_old[i][j] + temp_mat[i][j]/div;
    
    /*free up memory*/
    free(temp_arr);
    free_matrix(temp_mat,dim);
    
    return _Hi_new;
}

double ** math_tools::rep_reduce	(	char *	label,
		double **	rep_matrix,
		int	num_reps
	)			[protected]

Reduces a set of representation vectors.

Parameters:

	label	the point group label
	rep_matrix	the matrix of representations, each representation should be a row of this matrix whose length = number of irreps Parses input

Definition at line 24 of file symmetry.cc.

References psi::extrema::char_table::ctable, init_matrix(), psi::extrema::char_table::num_irreps, and psi::extrema::char_table::ops_coeffs.

                                                                            {
    
    int i,j, vec, rep, pg_order=0;
    char_table ctab(label);
    double **coef_matrix, red_coef;

    coef_matrix = init_matrix(num_reps,ctab.num_irreps);

    for(i=0;i<ctab.num_irreps;++i) 
        pg_order += ctab.ops_coeffs[i];

    int prod;
    for(rep=0;rep<num_reps;++rep) { 
        for(i=0;i<ctab.num_irreps;++i) {
            red_coef = 0.0;
            for(j=0;j<ctab.num_irreps;++j) {
                prod = ctab.ctable[i][j] * ctab.ops_coeffs[j];
                red_coef += rep_matrix[rep][j] * (double) prod;
            }
            coef_matrix[rep][i] = red_coef/(double)pg_order;
        }
    }

    return coef_matrix;
}

double ** math_tools::rep_project	(	char *	label,
		int	dim_vec,
		double **	result_vecs,
		int *	irrep_proj
	)			[protected]

Projects out irrep components from a reducible representation.

Parameters:

	label	point group label
	num_vars	dimension of the vector
	vectors	resulting from each symmetry operation
	irrep_proj	1 if irrep should be projected out and returned

Definition at line 62 of file symmetry.cc.

References psi::extrema::char_table::ctable, init_matrix(), psi::extrema::char_table::num_irreps, and psi::extrema::char_table::ops_coeffs.

                                                                        {
    
    int i,j,ir, proj_num, num_projections, pg_order;
    double **projected;
    char_table c(label);
    
    num_projections=0;
    for(ir=0;ir<c.num_irreps;++ir) 
        if( irrep_proj[ir] )
            ++num_projections;
    projected = init_matrix(num_projections,dim_vec);

    pg_order=0;
    for(i=0;i<c.num_irreps;++i) 
        pg_order += c.ops_coeffs[i];

    /* project out desired irreps */
    proj_num = 0;
    for(ir=0;ir<c.num_irreps;++ir) 
        if( irrep_proj[ir] ) {
            for(i=0;i<c.num_irreps;++i) 
                for(j=0;j<dim_vec;++j) {
                    projected[proj_num][j] += c.ctable[ir][i]*
                        result_vecs[i][j];
                }
            for(i=0;i<dim_vec;++i)
                projected[proj_num][i] /= (double) pg_order;
            ++proj_num;
        }

    return projected;
}

double ** math_tools::orthogonalize	(	int	nvecs,
		int	dimvecs,
		double **	vecs,
		int	normalize,
		double	norm_tol,
		int *	nindep
	)			[protected]

Othogonalizes a set of vectors.

Parameters:

	nvecs	number of vectors to orthogonalize
	dimvecs	length of each vector
	vec	nvec x dimvecs matrix containing vectors
	nindep	where number of independent vectors is saved
	norm_tol	norm of vectors big enough to keep
	normalize	1 if vectors should be normalized as well /*------------------------------------------------------------------------

Definition at line 111 of file symmetry.cc.

References ALMOST_ZERO, init_array(), and init_matrix().

                                                 {

    int i,j,k,p, num_big=0;
    double **omat_temp, **omat, norm, dot1, dot2, *temp;

    temp = init_array(dimvecs);

    /* max number of independent vectors is nvecs */
    omat_temp = init_matrix(nvecs,dimvecs);
    
    /* start with first vector */
    for(i=0;i<dimvecs;++i)
        omat_temp[0][i] = vecs[0][i]; 

    /* make each vector orthogonal to all previous vectors */
    for(i=1;i<nvecs;++i) {
        for(k=0;k<dimvecs;++k)
            omat_temp[i][k] = vecs[i][k]; 
        for(j=0;j<i;++j) {
            dot1=dot2=0.0;
            for(k=0;k<dimvecs;++k) {
                dot1 += vecs[i][k]*omat_temp[j][k];
                dot2 += omat_temp[j][k]*omat_temp[j][k];
            }
            if( (fabs(dot1)>ALMOST_ZERO) && (dot2>ALMOST_ZERO) ) {
                for(k=0;k<dimvecs;++k) 
                    omat_temp[i][k] -= dot1/dot2 * omat_temp[j][k];
            }
        }
    }

    /* if norm greater than tolerance, copy over to omat */
    for(i=0;i<nvecs;++i) {
        norm = 0.0;
        for(j=0;j<dimvecs;++j)
            norm += omat_temp[i][j]*omat_temp[i][j];
        if(norm>norm_tol) 
            ++num_big;
    }

    omat = init_matrix(num_big,dimvecs);

    p=0;
    for(i=0;i<nvecs;++i) {
        norm = 0.0;
        for(j=0;j<dimvecs;++j)
            norm += omat_temp[i][j]*omat_temp[i][j];
        if(norm>norm_tol) {
            for(j=0;j<dimvecs;++j)
                omat[p][j] = omat_temp[i][j];
            ++p;
        }
    }   

    (*nindep) = num_big;

    if(normalize) 
        for(i=0;i<num_big;++i) {
            norm=0;
            for(j=0;j<dimvecs;++j)
                norm += omat[i][j] * omat[i][j];
            for(j=0;j<dimvecs;++j)
                omat[i][j] /= sqrt(norm);
        }

    return omat;
}

---

The documentation for this class was generated from the following files:

• math_tools.h
• math_tools.cc
• symmetry.cc

---