libciomr: The PSI I/O and Math Library

Files
file	add_arr.cc
	Add two arrays.
file	add_mat.cc
	Add two matrices.
file	libciomr/block_matrix.cc
	Allocate a blocked (memory-contiguous) 2D matrix of doubles.
file	dot.cc
	Take dot product between two matrices.
file	eigout.cc
	Print eigenvectors and eigenvalues.
file	eigsort.cc
	Sort eigenvalues and eigenvectors in ascending or descending order.
file	eivout.cc
	Print eigenvectors and eigenvalues to output file.
file	ffile.cc
	Open PSI ASCII or small local binary (non-libpsio) files for reading/writing.
file	flin.cc
	Linear equation solver for A * x = b.
file	fndcor.cc
	Get the amount of core memory available from input.
file	init_matrix.cc
	Initialize a matrix of doubles.
file	int_array.cc
	This file includes the integer versions of several psi routines for handling arrays and matrices of doubles.
file	mmult.cc
	Multiply two matrices (superceded by C_DGEMM).
file	mxmb.cc
	Wrapper for the mmult function.
file	print_array.cc
	Print a lower-triangle array of doubles.
file	print_mat.cc
	Print a matrix of doubles.
file	psi_start.cc
	Initialize input, output, file prefix, etc.
file	psi_stop.cc
	Close input and output, stop input parser.
file	rsp.cc
	Diagonalize a symmetric matrix in packed (lower triangular) form.
file	sq_rsp.cc
	Diagnoalize a symmetrix square matrix.
file	sq_to_tri.cc
	Convert square matrix to lower triangle packing.
file	tqli.cc
	Diagonalizes a tridiagonal matrix output by tred2.
file	tred2.cc
	Converts a symmetric matrix to tridiagonal form for use in tqli.
file	tri_to_sq.cc
	Converts lower triangle to square matrix.
file	tstart.cc
	Controls starting and stopping of timers.
file	zero.cc
	Zero arrays or matrices of doubles.
Functions
void	add_arr (double *a, double *b, double *c, int n)
void	add_mat (double **a, double **b, double **c, int n, int m)
double **	block_matrix (unsigned long int n, unsigned long int m)
void	free_block (double **array)
double	dot_mat (double **a, double **b, int n)
void	dot_arr (double *a, double *b, int n, double *value)
void	eigout (double **a, double *b, double *c, int m, int n, FILE *out)
void	eigsort (double *d, double **v, int n)
void	mosort (double *d, double **v, int *sym, int nso, int nmo)
void	eivout (double **a, double *b, int m, int n, FILE *out)
void	ffile (FILE **fptr, char *suffix, int code)
void	ffile_noexit (FILE **fptr, char *suffix, int code)
void	ffileb (FILE **fptr, char *suffix, int code)
void	ffileb_noexit (FILE **fptr, char *suffix, int code)
void	flin (double **a, double *b, int in, int im, double *det)
void	fndcor (long int *maxcrb, FILE *infile, FILE *outfile)
double *	init_array (unsigned long int size)
double **	init_matrix (unsigned long int n, unsigned long int m)
void	free_matrix (double **array, unsigned long int size)
int *	init_int_array (int size)
void	zero_int_array (int *a, int size)
int **	init_int_matrix (int rows, int cols)
void	free_int_matrix (int **array)
void	zero_int_matrix (int **array, int rows, int cols)
void	print_int_mat (int **a, int m, int n, FILE *out)
void	mmult (double **AF, int ta, double **BF, int tb, double **CF, int tc, int nr, int nl, int nc, int add)
void	mxmb (double **a, int ia, int ja, double **b, int ib, int jb, double **c, int ic, int jc, int nrow, int nlnk, int ncol)
void	print_array (double *a, int m, FILE *out)
void	print_mat (double **a, int m, int n, FILE *out)
int	psi_start (FILE **infile, FILE **outfile, char **psi_file_prefix, int argc, char *argv[], int overwrite_output)
char *	psi_ifname ()
char *	psi_ofname ()
char *	psi_fprefix ()
int	psi_stop (FILE *infile, FILE *outfile, char *psi_file_prefix)
void	rsp (int nm, int n, int nv, double *array, double *e_vals, int matz, double **e_vecs, double toler)
void	sq_rsp (int nm, int n, double **array, double *e_vals, int matz, double **e_vecs, double toler)
void	sq_to_tri (double **bmat, double *amat, int size)
void	tqli (int n, double *d, double **z, double *e, int matz, double toler)
void	tred2 (int n, double **a, double *d, double *e, int matz)
void	tri_to_sq (double *amat, double **bmat, int size)
void	tstop (FILE *outfile)
void	zero_arr (double *a, int size)
void	zero_mat (double **a, int n, int m)

---

Detailed Description

---

Function Documentation

void add_arr	(	double *	a,
		double *	b,
		double *	c,
		int	n
	)

add_arr(): Add arrays a and b and put the result in array c. Adds the first n elements

Parameters:

	a	= first array to add
	b	= second array to add
	c	= array to hold the result of a+b
	n	= number of elements to add

Returns: none

Definition at line 23 of file add_arr.cc.

{
  register int i;

  for (i=0; i < n; i++) {
    c[i] = a[i]+b[i];
  }
}

void add_mat	(	double **	a,
		double **	b,
		double **	c,
		int	n,
		int	m
	)

add_mat(): Add matrices a and b into c for n rows and m columns

Parameters:

	a	= double star pointer to first matrix to add
	b	= double star pointer to second matrix to add
	c	= double star pointer to matrix to hold the result of a+b
	n	= number of rows in a,b,c
	m	= number of columns in a,b,c

Definition at line 20 of file add_mat.cc.

{
  register int i,j;

  if (n != m) {
    for (i=0; i < n ; i++) {
      for (j=0; j < m ; j++) {
        c[i][j] = a[i][j]+b[i][j];
      }
    }
  }
  else {
    for (i=0; i < n; i++) {
      for (j=0; j < i; j++) {
        c[i][j] = a[i][j]+b[i][j];
        c[j][i] = a[j][i]+b[j][i];
      }
      c[i][i] = a[i][i]+b[i][i];
    }
  }
}

double** block_matrix	(	unsigned long int	n,
		unsigned long int	m
	)

block_matrix(): Allocate a 2D array of doubles using contiguous memory

Allocates a contiguous block of memory for an array of doubles, allocates an array of pointers to the beginning of each row and returns the pointer to the first row pointer. This allows transparent 2d-array style access, but keeps memory together such that the matrix could be used in conjunction with FORTRAN matrix routines.

Allocates memory for an n x m matrix and returns a pointer to the first row.

Parameters:

	n	= number of rows (unsigned long to allow large matrices)
	m	= number of columns (unsigned long to allow large matrices)

Returns: double star pointer to newly allocated matrix

T. Daniel Crawford Sometime in 1994

Based on init_matrix() from libciomr

Definition at line 39 of file libciomr/block_matrix.cc.

Referenced by ael(), psi::detcas::bfgs_hessian(), psi::detcas::calc_gradient(), psi::detcas::calc_hessian(), david(), eri(), pople(), psi::extrema::zmat::write_chkpt(), and psi::extrema::coord_base_carts::write_chkpt().

   {
    double **A=NULL;
    double *B=NULL;
    unsigned long int i;

    if(!m || !n) return((double **) NULL);

    if ((A = (double **) malloc(n * (unsigned long int)sizeof(double *)))==NULL) {
         fprintf(stderr,"block_matrix: trouble allocating memory \n");
         fprintf(stderr,"n = %ld\n",n);
         exit(PSI_RETURN_FAILURE);
         }

    if ((B = (double *) malloc(m*n * (unsigned long int)sizeof(double)))==NULL) {
         fprintf(stderr,"block_matrix: trouble allocating memory \n");
         fprintf(stderr,"m = %ld\n",m);
         exit(PSI_RETURN_FAILURE);
         }

    bzero(B, m*n*(unsigned long int)sizeof(double));

    for (i = 0; i < n; i++) {
         A[i] = &(B[i*m]);
         }

    return(A);
   }

void dot_arr	(	double *	a,
		double *	b,
		int	n,
		double *	value
	)

dot_arr(): Take the dot product of the first n elements of two arrays a and b and return the result

Parameters:

	a	= first array to take dot product of
	b	= second array to take dot product of
	n	= number of elements in array
	value	= pointer to hold dot product result

Returns: none

Definition at line 51 of file dot.cc.

Referenced by normalize(), pople(), schmidt(), and schmidt_add().

{
  register int i;
  double tval;

  tval = 0.0;
  for (i=0; i < n; i++) {
    tval += a[i]*b[i];
  }
  *value = tval;
}

double dot_mat	(	double **	a,
		double **	b,
		int	n
	)

dot_mat(): Takes the dot product between two 2D matrices a and b with dimensions n x n and returns the value

Parameters:

	a	= first matrix for dot product
	b	= second matrix for dot product
	n	= number of rows/columns for matrices a and b

Returns: value of dot product

Definition at line 21 of file dot.cc.

{
  register int i,j;
  double *ta, *tb, tval;

  tval = 0.0;
  for (i=0; i < n; i++) {
    ta = a[i];
    tb = b[i];
    for (j=0; j < n; j++,ta++,tb++) {
      tval += (*ta) * (*tb);
    }
  }
  return(tval);
}

void eigout	(	double **	a,
		double *	b,
		double *	c,
		int	m,
		int	n,
		FILE *	out
	)

eigout(): Print out eigenvectors and eigenvalues.

Prints an n x m matrix of eigenvectors. Under each eigenvector, the corresponding elements of two arrays, b and c, will also be printed. This is useful for printing, for example, the SCF eigenvectors with their associated eigenvalues (orbital energies) and also the population.

Parameters:

	a	= matrix of eigenvectors (eigenvectors are columns)
	b	= first array to print under eigenvectors (e.g., eigenvalues)
	c	= second array to print under eigenvectors (e.g., populations)
	m	= number of rows in matrix a
	n	= number of columns in matrix a (and length of b and c)
	out	= file pointer for output

Returns: none

Definition at line 30 of file eigout.cc.

   {
      int ii,jj,kk,nn;
      int i,j;

      ii=0;jj=0;
L200:
      ii++;
      jj++;
      kk=10*jj;
      nn=n;
      if (nn > kk) nn=kk;
      fprintf (out,"\n");
      for (i=ii; i <= nn; i++) fprintf(out,"       %5d",i);
      fprintf (out,"\n");
      for (i=0; i < m; i++) {
         fprintf (out,"\n%5d",i+1);
         for (j=ii-1; j < nn; j++) {
            fprintf (out,"%12.7f",a[i][j]);
            }
         }
      fprintf (out,"\n");
      fprintf (out,"\n     ");
      for (j=ii-1; j < nn; j++) {
         fprintf(out,"%12.7f",b[j]);
         }
      fprintf (out,"\n");
      fprintf (out,"\n     ");
      for (j=ii-1; j < nn; j++) {
         fprintf(out,"%12.7f",c[j]);
         }
      fprintf (out,"\n");
      if (n <= kk) {
         fflush(out);
         return;
         }
      ii=kk; goto L200;
      }

void eigsort	(	double *	d,
		double **	v,
		int	n
	)

eigsort(): Sort the eigenvalues in d and eigenvectors in v in ascending (n>0) or descending (n<0) order. abs(n) is the number of eigenvalues.

Parameters:

	d	= array of eigenvalues
	v	= matrix of eigenvectors (each column is an eigenvector) Note: seems to assume v is a square matrix, could be a problem if, e.g., nmo != nso.
	n	= abs(n) is the number of eigenvalues/cols of v. Use n>0 to sort in ascending order, n<0 to sort in descending order

Returns: none

Definition at line 28 of file eigsort.cc.

Referenced by rsp(), and sq_rsp().

{
  int i,j,k;
  double p;

  /* Modified by Ed Valeev - if n is negative,
     sort eigenvalues in descending order */

  if (n >= 0) {
    for (i=0; i < n-1 ; i++) {
      k=i;
      p=d[i];
      for (j=i+1; j < n; j++) {
        if (d[j] < p) {
          k=j;
          p=d[j];
        }
      }
      if (k != i) {
        d[k]=d[i];
        d[i]=p;
        for (j=0; j < n; j++) {
          p=v[j][i];
          v[j][i]=v[j][k];
          v[j][k]=p;
        }
      }
    }
  }
  else {
    n = abs(n);
    for (i=0; i < n-1 ; i++) {
      k=i;
      p=d[i];
      for (j=i+1; j < n; j++) {
        if (d[j] > p) {
          k=j;
          p=d[j];
        }
      }
      if (k != i) {
        d[k]=d[i];
        d[i]=p;
        for (j=0; j < n; j++) {
          p=v[j][i];
          v[j][i]=v[j][k];
          v[j][k]=p;
        }
      }
    }
  }
}

void eivout	(	double **	a,
		double *	b,
		int	m,
		int	n,
		FILE *	out
	)

eivout: Print out eigenvectors and eigenvalues to the output file

Parameters:

	a	= eigenvectors
	b	= eigenvalues
	m	= rows of a
	n	= columns of a
	out	= output file pointer

Returns: none

Definition at line 24 of file eivout.cc.

{
  int ii,jj,kk,nn;
  int i,j;

  ii=0;jj=0;
L200:
  ii++;
  jj++;
  kk=10*jj;
  nn=n;
  if (nn > kk) nn=kk;
  fprintf (out,"\n");
  for (i=ii; i <= nn; i++) fprintf(out,"       %5d",i);
  fprintf (out,"\n");
  for (i=0; i < m; i++) {
    fprintf (out,"\n%5d",i+1);
    for (j=ii-1; j < nn; j++) {
      fprintf (out,"%12.7f",a[i][j]);
    }
  }
  fprintf (out,"\n");
  fprintf (out,"\n     ");
  for (j=ii-1; j < nn; j++) {
    fprintf(out,"%12.7f",b[j]);
  }
  fprintf (out,"\n");
  if (n <= kk) {
    fflush(out);
    return;
  }
  ii=kk; goto L200;
}

void ffile	(	FILE **	fptr,
		char *	suffix,
		int	code
	)

ffile(): Open a PSI3 ASCII file for reading/writing. Returns a pointer to the new file.

Parameters:

	suffix	= name of the file, not including automatic prefix
	code	= 0 (write), 1 (write/append), 2 (read)

Returns: none

Definition at line 28 of file ffile.cc.

Referenced by main(), psi::extrema::coord_base_carts::read_file11(), timer_done(), and psi::extrema::coord_base::write_opt().

{
  char name[100];

  /* build the standard file name */
  sprintf(name, "%s.%s", psi_file_prefix, suffix);

  switch (code) {
  case 0:
    *fptr = fopen(name,"w+");
    break;
  case 1:
    *fptr = fopen(name,"a+");
    break;
  case 2:
    *fptr = fopen(name,"r+");
    break;
  default:
    fprintf(stderr,"error in ffile: invalid code %d\n",code);
  }
  if (*fptr == NULL) { 
    fprintf(stderr,"error in ffile: cannot open file %s\n", suffix);
    exit(PSI_RETURN_FAILURE);
  }
}

void ffile_noexit	(	FILE **	fptr,
		char *	suffix,
		int	code
	)

ffile_noexit(): Open a PSI3 ASCII file for reading/writing. Returns a pointer to the new file via an argument. This function is the same as ffile(), but will not exit if fopen() fails.

Parameters:

	suffix	= name of the file, not including automatic prefix
	code	= 0 (write), 1 (write/append), 2 (read)

Returns: none

Definition at line 67 of file ffile.cc.

Referenced by psi::detcas::check_conv(), psi::extrema::deloc::deloc(), and psi::extrema::coord_base::read_opt().

{
  char name[100];

  /* build the standard file name */
  sprintf(name, "%s.%s", psi_file_prefix, suffix);

  switch (code) {
  case 0:
    *fptr = fopen(name,"w+");
    break;
  case 1:
    *fptr = fopen(name,"a+");
    break;
  case 2:
    *fptr = fopen(name,"r+");
    break;
  default:
    fprintf(stderr,"error in ffile_noexit: invalid code %d\n",code);
  }
}

void ffileb	(	FILE **	fptr,
		char *	suffix,
		int	code
	)

ffileb(): Open a PSI3 binary file for reading/writing. Returns a pointer to the new file.

Parameters:

	suffix	= name of the file, not including automatic prefix
	code	= 0 (write), 1 (write/append), 2 (read)

Returns: none

Definition at line 101 of file ffile.cc.

{
  char* name = (char*) malloc( (strlen(psi_file_prefix) + 
    strlen(suffix) + 2)*sizeof(char) );

  /* build the standard file name */
  sprintf(name, "%s.%s", psi_file_prefix, suffix);

  switch (code) {
  case 0:
    *fptr = fopen(name,"wb");
    break;
  case 1:
    *fptr = fopen(name,"ab");
    break;
  case 2:
    *fptr = fopen(name,"rb");
    break;
  default:
    fprintf(stderr,"error in ffileb: invalid code %d\n",code);
  }
  free(name);

  if (*fptr == NULL) {
    fprintf(stderr,"error in ffileb: cannot open file %s\n", suffix);
    exit(PSI_RETURN_FAILURE);
  }
}

void ffileb_noexit	(	FILE **	fptr,
		char *	suffix,
		int	code
	)

ffileb_noexit(): Open a PSI3 binary file for reading/writing. Returns a pointer to the new file via an argument. This function is the same as ffileb(), but will not exit if fopen() fails.

Parameters:

	suffix	= name of the file, not including automatic prefix
	code	= 0 (write), 1 (write/append), 2 (read)

Returns: none

Definition at line 143 of file ffile.cc.

{
  char* name = (char*) malloc( (strlen(psi_file_prefix) + 
    strlen(suffix) + 2)*sizeof(char) );

  /* build the standard file name */
  sprintf(name, "%s.%s", psi_file_prefix, suffix);

  switch (code) {
  case 0:
    *fptr = fopen(name,"wb");
    break;
  case 1:
    *fptr = fopen(name,"ab");
    break;
  case 2:
    *fptr = fopen(name,"rb");
    break;
  default:
    fprintf(stderr,"error in ffileb_noexit: invalid code %d\n",code);
  }
  free(name);
}

void flin	(	double **	a,
		double *	b,
		int	in,
		int	im,
		double *	det
	)

flin(): solves linear equations A * x = b.

Parameters:

	a	= coefficient matrix
	b	= known vectors
	in	= dimension of a(in*in)
	im	= number of b vectors
	det	= pointer to hold determinant of matrix a

Returns: none

Definition at line 29 of file flin.cc.

References init_array().

Referenced by pople().

{
  int i,j,k,*indx;

  indx = (int *) init_array(in);

  ludcmp(a,in,indx,det);

  for (i=0; i < in ; i++) *det *= a[i][i];

  for (j=0; j<im; j++)
    lubksb(a,in,indx,b+j*in);

  free(indx);
}

void fndcor	(	long int *	maxcrb,
		FILE *	infile,
		FILE *	outfile
	)

fndcor(): C translation of the Fortran version, to remove the need to link the library alloc, which also requires the linking of libparse, etc, etc...

This routine looks for the MEMORY keyword from input

Parameters:

	maxcrb	= long int ptr to hold size of maxcore in bytes
	infile	= file pointer to input file (eg input.dat)
	outfile	= file pointer to output file (eg output.dat)

Returns: none

David Sherrill, February 1994 Revised to handle more than 2GB of memory by Ed Valeev, October 2000

Revised to return the default if memory keyword is missing and raised the default to 256 MB. TDC, January 2003

Definition at line 43 of file fndcor.cc.

References ip_count(), ip_data(), ip_exist(), and ip_string().

Referenced by main().

{
  char type[20];
  char *s;
  int count;
  long int maxcrr;            /* maxcor in real words */
  char *maxcrr_str;           /* string representation of maxcrr */
  double size;
  int errcod; 

  maxcrr = DEF_MAXCRR;        /* set maxcor to default first */

  if(ip_exist("MEMORY",0)) { /* check if the keyword exists */
    errcod = ip_count("MEMORY", &count, 0);
    if (errcod != IPE_OK) fndcor_abort(infile, outfile);
    else if (errcod == IPE_NOT_AN_ARRAY) { /* Scalar specification of MEMORY */
      errcod = ip_string("MEMORY", &maxcrr_str, 0);
      if (errcod != IPE_OK) fndcor_abort(infile, outfile);
      maxcrr = atol(maxcrr_str);
    }
    /* Array specification of MEMORY */
    else if (count == 1) {
      errcod = ip_string("MEMORY", &maxcrr_str, 0);
      if (errcod != IPE_OK) fndcor_abort(infile, outfile) ;
      maxcrr = atol(maxcrr_str);
    }
    else if (count == 2) {
      errcod = ip_data("MEMORY", "%lf", &size, 1, 0);
      if (errcod != IPE_OK) fndcor_abort(infile, outfile);
      errcod = ip_data("MEMORY", "%s", type, 1, 1);
      if (errcod != IPE_OK) fndcor_abort(infile, outfile);
      /* convert string to uppercase */
      for (s=type; *s!='\0'; s++) {
        if (*s>='a' && *s <='z') *s = *s + 'A' - 'a';
      }
      if ((strcmp(type, "R")==0) || (strcmp(type, "REAL")==0))
        maxcrr = (long int) size;
      else if ((strcmp(type, "I")==0) || (strcmp(type, "INTEGER")==0)) 
        maxcrr = (long int) (size * sizeof(int) / sizeof(double));
      else if ((strcmp(type, "B")==0) || (strcmp(type, "BYTES")==0)) 
        maxcrr = (long int) (size / sizeof(double));
      else if ((strcmp(type, "KB")==0) || (strcmp(type, "KBYTES")==0)) 
        maxcrr = (long int) (1000.0 * size / sizeof(double));
      else if ((strcmp(type, "MB")==0) || (strcmp(type, "MBYTES")==0))
        maxcrr = (long int) (1000000.0 * size / sizeof(double));
      else if ((strcmp(type, "GB")==0) || (strcmp(type, "GBYTES")==0))
        maxcrr = (long int) (1000000000.0 * size / sizeof(double));
      else {
        fprintf(outfile, "bad data type, specify one of: \n") ;
        fprintf(outfile, "REAL, INTEGER, BYTES, KBYTES, MBYTES, or GBYTES\n");
        fndcor_abort(infile, outfile);
      }
    }
  }
  
  *maxcrb = maxcrr * sizeof(double);

  return;
}

void free_block

(

double **

array

)

free_block(): Free a block matrix

Parameters:

array

= pointer to matrix to be freed

Returns: none

Definition at line 78 of file libciomr/block_matrix.cc.

   {
     if(array == NULL) return;
      free(array[0]);
      free(array);
   }

void free_int_matrix

(

int **

array

)

free_int_matrix(): Free a matrix of integers. Pass a pointer to the matrix.

Parameters:

array

= pointer to integer matrix

Definition at line 111 of file int_array.cc.

Referenced by ras_set(), ras_set2(), and psi::extrema::deloc::~deloc().

{
  free(array[0]);
  free(array);
}

void free_matrix	(	double **	array,
		unsigned long int	size
	)

free_matrix(): Free a 2D matrix allocated with init_matrix().

Parameters:

	array	= matrix to free
	size	= number of rows (unsigned long to allow large matrices)

Returns: none

Definition at line 61 of file init_matrix.cc.

Referenced by psi::cscf::check_rot(), psi::extrema::internals::compute_A(), psi::extrema::internals::compute_G(), psi::extrema::deloc::deloc(), psi::extrema::internals::grad_trans(), psi::extrema::coord_base::H_test(), mmult(), psi::extrema::zmat::print_c_grads(), psi::extrema::coord_base_carts::print_c_grads(), psi::extrema::zmat::print_carts(), psi::extrema::coord_base_carts::print_carts(), psi::extrema::coord_base::print_H(), sq_rsp(), psi::extrema::math_tools::update_bfgs(), psi::extrema::math_tools::update_ms(), psi::extrema::zmat::zmat(), psi::extrema::coord_base::~coord_base(), and psi::extrema::internals::~internals().

{
  unsigned long int i;

  for (i=0; i < size ; i++) {
    free(array[i]);
  }

  free(array);
}

double* init_array

(

unsigned long int

size

)

init_array(): This function initializes an array of doubles of length 'size' and returns a pointer to the first element

Parameters:

size

= length of array (unsigned long to allow large arrays)

Returns: pointer to new array

Definition at line 24 of file init_array.cc.

Referenced by ael(), psi::extrema::internals::B_row_angle(), psi::extrema::internals::B_row_bond(), psi::extrema::internals::B_row_tors(), psi::extrema::zmat::back_transform(), psi::extrema::internals::back_transform(), psi::detcas::bfgs_hessian(), psi::detcas::calc_gradient(), psi::detcas::calc_hessian(), psi::extrema::coord_base_carts::coord_base_carts(), david(), psi::extrema::deloc::deloc(), psi::detcas::ds_hessian(), eri(), flin(), psi::extrema::coord_base::H_test(), invert_matrix(), psi::extrema::internals::mem_alloc(), psi::extrema::coord_base::mem_alloc(), psi::extrema::zmat::newton_step(), psi::extrema::math_tools::newton_step(), norm_const(), psi::extrema::math_tools::orthogonalize(), pople(), psi::detcas::rotate_orbs(), rsp(), psi::detcas::scale_gradient(), schmidt(), schmidt_add(), slaterdetvector_init(), sq_rsp(), psi::extrema::coord_base::update_Hi(), psi::extrema::math_tools::update_ms(), and psi::extrema::zmat::zmat().

{
  double *array;

  if ((array = (double *) malloc(size*(unsigned long int)sizeof(double)))
    == NULL) {
    fprintf(stderr,"init_array: trouble allocating memory \n");
    fprintf(stderr,"size = %ld\n",size);
    exit(PSI_RETURN_FAILURE);
  }
  bzero(array,size*(unsigned long int)sizeof(double));
  return(array);
}

int* init_int_array

(

int

size

)

init_int_array(): Allocates memory for one-D array of ints of dimension 'size' and returns pointer to 1st element. Zeroes all elements.

Just modified the init_array() routine to do int's instead. This will avoid the temptation to allocate 5 integers by p = (int *) init_array(5/2), which is bad.

Parameters:

size

= length of array to allocate

Returns: pointer to new array

C. David Sherrill

Definition at line 34 of file int_array.cc.

Referenced by psi::detcas::bfgs_hessian(), david(), psi::extrema::deloc::deloc(), get_frzcpi(), get_frzvpi(), invert_matrix(), main(), ras_set(), ras_set2(), reorder_qt(), reorder_qt_uhf(), and psi::extrema::zmat::zmat().

{
  int *array;

  if ((array = (int *) malloc(sizeof(int)*size))==NULL) {
    fprintf(stderr,"init_array:  trouble allocating memory \n");
    fprintf(stderr,"size = %d\n",size);
    exit(PSI_RETURN_FAILURE);
  }
  bzero(array,sizeof(int)*size);
  return(array);
}

int** init_int_matrix	(	int	rows,
		int	cols
	)

init_int_matrix(): Function initializes (allocates and clears) a matrix of integers with dimensions 'rows' by 'cols' and returns a pointer to it (ptr to first row ptr). The matrix layout is blocked, i.e. like produced by block_matrix()

Parameters:

	rows	= number of rows
	cols	= number of columns

Returns: pointer to first row of newly-allocated integer block matrix

Definition at line 78 of file int_array.cc.

Referenced by ras_set(), and ras_set2().

{
   int **array=NULL;
   int i;

   if ((array = (int **) malloc(sizeof(int *)*rows))==NULL) {
     fprintf(stderr,"init_int_matrix: trouble allocating memory \n"); 
     fprintf(stderr,"rows = %d\n", rows);
     exit(PSI_RETURN_FAILURE);
   }

   if ((array[0] = (int *) malloc (sizeof(int)*cols*rows))==NULL) {
     fprintf(stderr,"init_int_matrix: trouble allocating memory \n"); 
     fprintf(stderr,"rows = %d, cols = %d", rows, cols);
     exit(PSI_RETURN_FAILURE) ;
   }
   for (i=1; i<rows; i++) {
     array[i] = array[i-1] + cols;
   }
   bzero(array[0], sizeof(int)*cols*rows);

   return array;
}

double** init_matrix	(	unsigned long int	n,
		unsigned long int	m
	)

init_matrix(): Initialize an nxm matrix of doubles and return a pointer to the first row. Note that this does not form a matrix which is necessarily contiguous in memory. Use block_matrix() for that.

Parameters:

	n	= number of rows (unsigned long to allow large matrices)
	m	= number of columns (unsigned long to allow large matrices)

Returns: pointer to first row

Definition at line 26 of file init_matrix.cc.

Referenced by psi::cscf::check_rot(), psi::extrema::internals::compute_A(), psi::extrema::internals::compute_G(), psi::extrema::deloc::deloc(), psi::extrema::internals::grad_trans(), psi::extrema::coord_base::H_test(), psi::extrema::internals::mem_alloc(), psi::extrema::coord_base::mem_alloc(), mmult(), psi::extrema::math_tools::orthogonalize(), psi::extrema::zmat::print_c_grads(), psi::extrema::coord_base_carts::print_c_grads(), psi::extrema::zmat::print_carts(), psi::extrema::coord_base_carts::print_carts(), psi::extrema::math_tools::rep_project(), psi::extrema::math_tools::rep_reduce(), sq_rsp(), psi::extrema::math_tools::update_bfgs(), and psi::extrema::math_tools::update_ms().

{
  double **array=NULL;
  unsigned long int i;

  if ((array = (double **) malloc(n*(unsigned long int)sizeof(double *)))
    ==NULL) {
    fprintf(stderr,"init_matrix: trouble allocating memory \n");
    fprintf(stderr,"n = %ld\n",n);
    exit(PSI_RETURN_FAILURE);
  }

  for (i = 0; i < n; i++) {
    if ((array[i] = (double *) malloc(m*(unsigned long int)sizeof(double)))
      ==NULL) {
      fprintf(stderr,"init_matrix: trouble allocating memory \n");
      fprintf(stderr,"i = %ld m = %ld\n",i,m);
      exit(PSI_RETURN_FAILURE);
    }
    bzero(array[i],m*(unsigned long int)sizeof(double));
  }
  return(array);
}

void mmult	(	double **	AF,
		int	ta,
		double **	BF,
		int	tb,
		double **	CF,
		int	tc,
		int	nr,
		int	nl,
		int	nc,
		int	add
	)

mmult(): a reasonably fast matrix multiply (at least on the DEC3100) written by ETS

Parameters:

	AF	= first matrix to multiply
	ta	= if 1, transpose AF before multiplying; otherwise, 0
	BF	= second matrix to multiply
	tb	= if 1, transpose BF before multiplying; otherwise 0
	CF	= matrix to hold result of AF*BF
	tc	= if 1, transpose CF after the multiplication; otherwise 0
	nr	= number of rows of AF
	nl	= number of cols of AF and rows of BF
	nc	= number of cols of BF
	add	= if 1, add AF*BF to the matrix passed in as CF; else 0

Returns: none

nr,nl,nc are the number of rows,links,and columns in the final matrices to be multiplied together if ta=0 AF should have the dimensions nr x nl if ta=1 AF should have the dimensions nl x nr if tb=0 BF should have the dimensions nl x nc if tb=1 BF should have the dimensions nc x nl if tc=0 CF should have the dimensions nr x nc if tc=1 CF should have the dimensions nc x nr

Definition at line 46 of file mmult.cc.

References free_matrix(), and init_matrix().

Referenced by psi::cscf::check_rot(), psi::extrema::internals::compute_A(), psi::extrema::internals::compute_G(), psi::extrema::deloc::deloc(), psi::extrema::internals::grad_trans(), mxmb(), and psi::extrema::math_tools::update_bfgs().

{
   int odd_nr,odd_nc,odd_nl;
   int i,j,k,ij;
   double t00,t01,t10,t11;
   double **a,**b;
   double *att,*bt;
   double *at1,*bt1;

   if(!aa) {
      aa = (double **) init_matrix(nr,nl);
      bb = (double **) init_matrix(nc,nl);
      keep_nr = nr;
      keep_nl = nl;
      keep_nc = nc;
      }

   if(nl > keep_nl) {
      free_matrix(aa,keep_nr);
      free_matrix(bb,keep_nc);
      keep_nl = nl;
      keep_nr = (nr > keep_nr) ? nr : keep_nr;
      keep_nc = (nc > keep_nc) ? nc : keep_nc;
      aa = (double **) init_matrix(keep_nr,keep_nl);
      bb = (double **) init_matrix(keep_nc,keep_nl);
      }
   if(nr > keep_nr) {
      free_matrix(aa,keep_nr);
      keep_nr = nr;
      aa = (double **) init_matrix(keep_nr,keep_nl);
      }
   if(nc > keep_nc) {
      free_matrix(bb,keep_nc);
      keep_nc = nc;
      bb = (double **) init_matrix(keep_nc,keep_nl);
      }

   odd_nr = (nr)%2;
   odd_nc = (nc)%2;
   odd_nl = (nl)%2;

   a=aa;
   if(ta)
      for(i=0; i < nr ; i++)
         for(j=0; j < nl ; j++)
            a[i][j] = AF[j][i];
   else
      a=AF;

   b=bb;
   if(tb)
      b=BF;
   else
      for(i=0; i < nc ; i++)
         for(j=0; j < nl ; j++)
            b[i][j] = BF[j][i];
      
   for(j=0; j < nc-1 ; j+=2) {
      for(i=0; i < nr-1 ; i+=2) {
         att=a[i]; bt=b[j];
         at1=a[i+1]; bt1=b[j+1];
         if(add) {
            if(tc) {
               t00 = CF[j][i];
               t01 = CF[j+1][i];
               t10 = CF[j][i+1];
               t11 = CF[j+1][i+1];
               }
            else {
               t00 = CF[i][j];
               t01 = CF[i][j+1];
               t10 = CF[i+1][j];
               t11 = CF[i+1][j+1];
               }
            }
         else t00=t01=t10=t11=0.0;
         for(k=nl; k ; k--,att++,bt++,at1++,bt1++) {
            t00 += *att * *bt;
            t01 += *att * *bt1;
            t10 += *at1 * *bt;
            t11 += *at1 * *bt1;
            }
         if(tc) {
            CF[j][i]=t00;
            CF[j+1][i]=t01;
            CF[j][i+1]=t10;
            CF[j+1][i+1]=t11;
            }
         else {
            CF[i][j]=t00;
            CF[i][j+1]=t01;
            CF[i+1][j]=t10;
            CF[i+1][j+1]=t11;
            }
         }
      if(odd_nr) {
         att=a[i]; bt=b[j];
         bt1=b[j+1];
         if(add) {
            if(tc) {
               t00 = CF[j][i];
               t01 = CF[j+1][i];
               }
            else {
               t00 = CF[i][j];
               t01 = CF[i][j+1];
               }
            }
         else t00=t01=0.0;
         for(k= nl; k ; k--,att++,bt++,bt1++) {
            t00 += *att * *bt;
            t01 += *att * *bt1;
            }
         if(tc) {
            CF[j][i]=t00;
            CF[j+1][i]=t01;
            }
         else {
            CF[i][j]=t00;
            CF[i][j+1]=t01;
            }
         }
      }
   if(odd_nc) {
      for(i=0; i < nr-1 ; i+=2) {
         att=a[i]; bt=b[j];
         at1=a[i+1];
         if(add) {
            if(tc) {
               t00 = CF[j][i];
               t10 = CF[j][i+1];
               }
            else {
               t00 = CF[i][j];
               t10 = CF[i+1][j];
               }
            }
         else t00=t10=0.0;
         for(k= nl; k ; k--,att++,bt++,at1++) {
            t00 += *att * *bt;
            t10 += *at1 * *bt;
            }
         if(tc) {
            CF[j][i]=t00;
            CF[j][i+1]=t10;
            }
         else {
            CF[i][j]=t00;
            CF[i+1][j]=t10;
            }
         }
      if(odd_nr) {
         att=a[i]; bt=b[j];
         if(add)
            t00 = (tc) ? CF[j][i] : CF[i][j];
         else t00=0.0;
         for(k=nl; k ; k--,att++,bt++)
            t00 += *att * *bt;
         if(tc) CF[j][i]=t00;
         else CF[i][j]=t00;
         }
      }
   }

void mosort	(	double *	d,
		double **	v,
		int *	sym,
		int	nso,
		int	nmo
	)

mosort(): Minor modification of eigsort() to also sort a series of irrep labels.

Parameters:

	d	= array of eigenvalues
	v	= matrix of eigenvectors (each column is an eigenvector)
	sym	= array of symmetry ID's (irreps)
	nso	= number of rows in v
	nmo	= abs(nmo) is the number of eigenvalues/cols of v. Use nmo>0 to sort in ascending order, nmo<0 to sort in descending order

Returns:none

TDC, 6/03

Definition at line 99 of file eigsort.cc.

{
  int i, j, k, l;
  double p;

  if(nmo > 0) {
    for (i=0; i < nmo-1 ; i++) {
      k=i;
      p=d[i];
      for (j=i+1; j < nmo; j++) {
        if (d[j] < p) {
          k=j;
          p=d[j];
        }
      }
      if (k != i) {
        d[k]=d[i];
        d[i]=p;

        l = sym[i];
        sym[i] = sym[k];
        sym[k] = l;

        for (j=0; j < nso; j++) {
          p=v[j][i];
          v[j][i]=v[j][k];
          v[j][k]=p;
        }
      }
    }
  }
  else if(nmo < 0) {
    nmo = abs(nmo);
    for (i=0; i < nmo-1 ; i++) {
      k=i;
      p=d[i];
      for (j=i+1; j < nmo; j++) {
        if (d[j] > p) {
          k=j;
          p=d[j];
        }
      }
      if (k != i) {
        d[k]=d[i];
        d[i]=p;

        l = sym[i];
        sym[i] = sym[k];
        sym[k] = l;

        for (j=0; j < nso; j++) {
          p=v[j][i];
          v[j][i]=v[j][k];
          v[j][k]=p;
        }
      }
    }
  }
}

void mxmb	(	double **	a,
		int	ia,
		int	ja,
		double **	b,
		int	ib,
		int	jb,
		double **	c,
		int	ic,
		int	jc,
		int	nrow,
		int	nlnk,
		int	ncol
	)

mxmb: multiplies two rectangular matrices together (wrapper for mmult). Deprecated; use C_DGEMM instead.

Parameters:

	a	= first matrix to multiply
	ia	= if 1, normal multiplication of a
	ja	= if 1, transpose a before multiplication
	b	= second matrix to multiply
	ib	= if 1, normal multiplication of b
	jb	= if 1, transpose b before multiplication
	c	= matrix to store the result
	ic	= if 1, normal multiplication into c
	jb	= if 1, transpose c after multiplication
	nrow	= number of rows of a
	nlnk	= number of columns of a and rows of b
	ncol	= number of columns of b

Returns: none

Definition at line 41 of file mxmb.cc.

References mmult().

   {
      if (ic == 1) {
         if (ia == 1) {
            if (ib == 1) {
               mmult(a,0,b,0,c,0,nrow,nlnk,ncol,0);
               }
            else {
               if (jb == 1) {
                  mmult(a,0,b,1,c,0,nrow,nlnk,ncol,0);
                  }
               else {
                  mxmbol(a,ia,ja,b,ib,jb,c,ic,jc,nrow,nlnk,ncol);
                  }
               }
            }
         else {
            if (ja == 1) {
               if (ib == 1) {
                  mmult(a,1,b,0,c,0,nrow,nlnk,ncol,0);
                  }
               else {
                  if (jb == 1) {
                     mmult(a,1,b,1,c,0,nrow,nlnk,ncol,0);
                     }
                  else {
                     mxmbol(a,ia,ja,b,ib,jb,c,ic,jc,nrow,nlnk,ncol);
                     }
                  }
               }
            else {
               mxmbol(a,ia,ja,b,ib,jb,c,ic,jc,nrow,nlnk,ncol);
               }
            }
         }
      else {
         if (jc == 1) {
            if (ia == 1) {
               if (ib == 1) {
                  mmult(a,0,b,0,c,1,nrow,nlnk,ncol,0);
                  }
               else {
                  if (jb == 1) {
                     mmult(a,0,b,1,c,1,nrow,nlnk,ncol,0);
                     }
                  else {
                     mxmbol(a,ia,ja,b,ib,jb,c,ic,jc,nrow,nlnk,ncol);
                     }
                  }
               }
            else {
               if (ja == 1) {
                  if (ib == 1) {
                     mmult(a,1,b,0,c,1,nrow,nlnk,ncol,0);
                     }
                  else {
                     if (jb == 1) {
                        mmult(a,1,b,1,c,1,nrow,nlnk,ncol,0);
                        }
                     else {
                        mxmbol(a,ia,ja,b,ib,jb,c,ic,jc,nrow,nlnk,ncol);
                        }
                     }
                  }
               else {
                  mxmbol(a,ia,ja,b,ib,jb,c,ic,jc,nrow,nlnk,ncol);
                  }
               }
            }
         else {
            mxmbol(a,ia,ja,b,ib,jb,c,ic,jc,nrow,nlnk,ncol);
            }
         }
      }

void print_array	(	double *	a,
		int	m,
		FILE *	out
	)

print_array(): Prints a lower-triangle of a symmetric matrix packed as an array of doubles.

Parameters:

	a	= array (packed lower triangle of matrix) to print
	m	= dimension of matrix (mxm)
	out	= file pointer for output

Returns: none

Definition at line 23 of file print_array.cc.

Referenced by iwl_rdone().

   {
      int ii,jj,kk,mm,nn,ll;
      int i,j,k,i1,i2;

      ii=0;jj=0;
L200:
      ii++;
      jj++;
      kk=10*jj;
      nn = kk + kk*(kk-1)/2;
      mm=m;
      if (m > kk) mm=kk;
      ll = 2*(mm-ii+1)+1;
      fprintf (out,"\n");
      for (i=ii; i <= mm; i++) fprintf(out,"       %5d",i);
      fprintf (out,"\n");
      for (i=ii; i <= m; i++) {
         i1=i*(i-1)/2+ii;
         i2=i+i*(i-1)/2;
         if (i2 > nn) i2 = i1+9;
         fprintf (out,"\n%5d",i);
         for (j=i1; j <= i2; j++) {
            fprintf (out,"%12.7f",a[j-1]);
            }
         }
      if (m <= kk) {
         fprintf(out,"\n");
         fflush(out);
         return;
         }
      ii=kk; goto L200;
      }

void print_int_mat	(	int **	a,
		int	m,
		int	n,
		FILE *	out
	)

print_int_mat(): Print a matrix of integers. Pass the matrix, the number of rows and columns, and the output file pointer.

Parameters:

	a	= integer matrix to print
	m	= number of rows in matrix
	n	= number of columns in matrix
	out	= FILE pointer to output file

Returns: none

Definition at line 151 of file int_array.cc.

{
  int ii,jj,kk,nn,ll;
  int i,j,k;

  ii=0;jj=0;
L200:
  ii++;
  jj++;
  kk=10*jj;
  nn=n;
  if (nn > kk) nn=kk;
  ll = 2*(nn-ii+1)+1;
  fprintf (out,"\n   ");
  for (i=ii; i <= nn; i++) fprintf(out,"   %5d",i);
  fprintf (out,"\n");
  for (i=0; i < m; i++) {
    fprintf (out,"\n%5d",i+1);
    for (j=ii-1; j < nn; j++) {
      fprintf (out,"%8d",a[i][j]);
    }
  }
  fprintf (out,"\n");
  if (n <= kk) {
    fflush(out);
    return;
  }
  ii=kk; goto L200;
}

void print_mat	(	double **	a,
		int	m,
		int	n,
		FILE *	out
	)

print_mat: Print a matrix a of dimensions mxn to file pointer out.

Parameters:

	a	= matrix to print
	m	= number of rows in matrix
	n	= number of columns in matrix
	out	= file pointer for output

Returns: none

Definition at line 23 of file print_mat.cc.

   {
      int ii,jj,kk,nn,ll;
      int i,j,k;

      ii=0;jj=0;
L200:
      ii++;
      jj++;
      kk=10*jj;
      nn=n;
      if (nn > kk) nn=kk;
      ll = 2*(nn-ii+1)+1;
      fprintf (out,"\n");
      for (i=ii; i <= nn; i++) fprintf(out,"       %5d",i);
      fprintf (out,"\n");
      for (i=0; i < m; i++) {
         fprintf (out,"\n%5d",i+1);
         for (j=ii-1; j < nn; j++) {
            fprintf (out,"%12.7f",a[i][j]);
            }
         }
      fprintf (out,"\n");
      if (n <= kk) {
         fflush(out);
         return;
         }
      ii=kk; goto L200;
      }

char* psi_fprefix

(

)

psi_fprefix()

This function returns the PSI file prefix

Arguments: none

Returns: the pointer to the string containing the PSI file prefix if it has been determined, NULL otherwise

Definition at line 291 of file psi_start.cc.

{
  return fprefix;
}

char* psi_ifname

(

)

psi_ifname()

This function returns the input file name

Arguments: none

Returns: the pointer to the string containing the input file name if it has been determined, NULL otherwise

Definition at line 257 of file psi_start.cc.

{
  return ifname;
}

char* psi_ofname

(

)

psi_ofname()

This function returns the output file name

Arguments: none

Returns: the pointer to the string containing the output file name if it has been determined, NULL otherwise

Definition at line 274 of file psi_start.cc.

{
  return ofname;
}

int psi_start	(	FILE **	infile,
		FILE **	outfile,
		char **	psi_file_prefix,
		int	argc,
		char *	argv[],
		int	overwrite_output
	)

psi_start(): This function initializes the input, output files, file prefix, etc., by checking command line arguments and environmental variables. It also initializes the Input Parsing library.

Parameters:

	argc	= number of command-line arguments passed
	argv	= command-line arguments
	overwrite	= whether to overwrite output file (1) or append to it (0). Most PSI modules will want to append.

Returns: one of standard PSI error codes

Definition at line 38 of file psi_start.cc.

References ip_string().

{
  int i, errcod;
                                /* state flags */
  int found_if_np = 0;          /* found input file name without -f */
  int found_of_np = 0;          /* found output file name without -o */
  int found_fp_np = 0;          /* found file prefix name without -p */
  int found_if_p = 0;           /* found input file name with -f */
  int found_of_p = 0;           /* found output file name with -o */
  int found_fp_p = 0;           /* found file prefix name with -p */
  char *cfname = NULL;
  char *userhome;
  FILE *psirc;
  char *arg;
  char *tmpstr1;

  /* process command-line arguments in sequence */
  for(i=0; i<argc; i++) {
    arg = argv[i];
    if (!strcmp(arg,"-f") && !found_if_p) {
      ifname = argv[++i];
      found_if_p = 1;
    }
    else if (!strcmp(arg,"-o") && !found_of_p) {
      ofname = argv[++i];
      found_of_p = 1;
    }
    else if (!strcmp(arg,"-p") && !found_fp_p) {
      fprefix = argv[++i];
      found_fp_p = 1;
    }
    else if (arg[0] == '-') {
      fprintf(stderr, "Error: unrecognized command-line argument %s\n", arg);
      return(PSI_RETURN_FAILURE);
    }
    else if (!found_if_np) {
      ifname = arg;
      found_if_np = 1;
    }
    else if (!found_of_np) {
      ofname = arg;
      found_of_np = 1;
    }
    else if (!found_fp_np) {
      fprefix = arg;
      found_fp_np = 1;
    }
    else {
      fprintf(stderr, "Error: too many command-line arguments given\n");
      return(PSI_RETURN_FAILURE);
    }
  }
  

  /* check if some args were specified in both prefixed and nonprefixed form */
  if (found_if_p && found_if_np) {
    fprintf(stderr, 
      "Error: input file name specified both with and without -f\n");
    fprintf(stderr, 
      "Usage: (module) [options] -f input -o output [-p prefix]  OR\n");
    fprintf(stderr, "       (module) [options] input output [prefix]\n");
    return(PSI_RETURN_FAILURE);
  }
  if (found_of_p && found_of_np) {
    fprintf(stderr, 
      "Error: output file name specified both with and without -o\n");
    fprintf(stderr, 
      "Usage: (module) [options] -f input -o output [-p prefix]  OR\n");
    fprintf(stderr, "       (module) [options] input output [prefix]\n");
    return(PSI_RETURN_FAILURE);
  }
  if (found_fp_p && found_fp_np) {
    fprintf(stderr, 
      "Error: file prefix specified both with and without -p\n");
    fprintf(stderr, 
      "Usage: (module) [options] -f input -o output -p prefix  OR\n");
    fprintf(stderr, "       (module) [options] input output prefix\n");
    return(PSI_RETURN_FAILURE);
  }

  
  /* if some args were not specified on command-line - check the environment */
  if (ifname == NULL)
    ifname = getenv("PSI_INPUT");
  if (ofname == NULL)
    ofname = getenv("PSI_OUTPUT");
  if (fprefix == NULL)
    fprefix = getenv("PSI_PREFIX");

  /* if some arguments still not defined - assign default values */
  if (ifname == NULL)
    ifname = strdup("input.dat");
  if (ofname == NULL)
    ofname = strdup("output.dat");
  /* default prefix is not assigned here yet because need to check 
     input file first */

  /* open input and output files */
  if(ifname[0]=='-' && ifname[1]=='\x0')
    *infile=stdin;
  else
    *infile = fopen(ifname, "r");
  if (*infile == NULL) {
    fprintf(stderr, "Error: could not open input file %s\n",ifname);
    return(PSI_RETURN_FAILURE);
  }
  if (overwrite_output)
  {
    if(ofname[0]=='-' && ofname[1]=='\x0')
      *outfile=stdout;
    else
      *outfile = fopen(ofname, "w");
  }
  else
  {
    if(ofname[0]=='-' && ofname[1]=='\x0')
      *outfile=stdout;
    else
      *outfile = fopen(ofname, "a");
  }
  if (*outfile == NULL) {
    fprintf(stderr, "Error: could not open output file %s\n",ofname);
    return(PSI_RETURN_FAILURE);
  }

  /* initialize libipv1 */
  ip_set_uppercase(1);
  ip_initialize(*infile, *outfile);
  ip_cwk_clear();

  /* open user's PSI configuration file (default, $HOME/.psirc) */
  cfname = getenv("PSI_RC");
  if (cfname == NULL) {
    userhome = getenv("HOME");
    cfname = (char *) malloc((10+strlen(userhome))*sizeof(char));
    sprintf(cfname, "%s%s", userhome, "/.psirc");
    psirc = fopen(cfname, "r");
    free(cfname);
  }
  else psirc = fopen(cfname, "r");
  if(psirc != NULL) {
    ip_append(psirc, stderr);
    fclose(psirc);
  }  

  /* lastly, everybody needs DEFAULT and PSI sections */
  ip_cwk_add(":DEFAULT");
  ip_cwk_add(":PSI");

  /* if prefix still NULL - check input file */
  if (fprefix == NULL)
    errcod = ip_string(":DEFAULT:FILES:DEFAULT:NAME",&fprefix,0);
  if (fprefix == NULL)
    errcod = ip_string(":DEFAULT:NAME",&fprefix,0);
  if (fprefix == NULL)
    errcod = ip_string(":PSI:FILES:DEFAULT:NAME",&fprefix,0);
  if (fprefix == NULL)
    errcod = ip_string(":PSI:NAME",&fprefix,0);

  /* copy over file prefix, etc. into their appropriate variables */
  if (fprefix == NULL) {
    fprefix = strdup(PSI_DEFAULT_FILE_PREFIX);
  }
  *psi_file_prefix = strdup(fprefix);

  /* other Psi modules called by this module should read from the same input 
     file set the value of PSI_INPUT for the duration of this run */
#if HAVE_PUTENV
  tmpstr1 = (char *) malloc(11+strlen(ifname));
  sprintf(tmpstr1, "PSI_INPUT=%s", ifname);
  putenv(tmpstr1);  /* note potential memory leak */
#elif HAVE_SETENV
  setenv("PSI_OUTPUT",ifname,1);
#else
#error "Have neither putenv nor setenv. Something must be very broken on this system."
#endif

  /* By default, other Psi modules called by this module should write to 
     the same output file set the value of PSI_OUTPUT for the duration of 
     this run */
#if HAVE_PUTENV
  tmpstr1 = (char *) malloc(12+strlen(ofname));
  sprintf(tmpstr1, "PSI_OUTPUT=%s", ofname);
  putenv(tmpstr1); /* note potential memory leak */
#elif HAVE_SETENV
  setenv("PSI_OUTPUT",ofname,1);
#else
#error "Have neither putenv nor setenv. Something must be very broken on this system."
#endif

  /* By default, other Psi modules called by this module should use the same
     prefix too set the value of PSI_PREFIX for the duration of this run */
#if HAVE_PUTENV
  tmpstr1 = (char *) malloc(12+strlen(fprefix));
  sprintf(tmpstr1, "PSI_PREFIX=%s", fprefix);
  putenv(tmpstr1); /* note potential memory leak */
#elif HAVE_SETENV
  setenv("PSI_PREFIX",fprefix,1);
#else
#error "Have neither putenv nor setenv. Something must be very broken on this system."
#endif

  return(PSI_RETURN_SUCCESS);
}

int psi_stop	(	FILE *	infile,
		FILE *	outfile,
		char *	psi_file_prefix
	)

psi_stop(): This function closes input and output files and deinitializes Input Parsing library.

Arguments: none

Returns: one of standard PSI error codes

Definition at line 26 of file psi_stop.cc.

{
  ip_done();
  free(psi_file_prefix);
  fclose(outfile);
  fclose(infile);

  return(PSI_RETURN_SUCCESS);
}

void rsp	(	int	nm,
		int	n,
		int	nv,
		double *	array,
		double *	e_vals,
		int	matz,
		double **	e_vecs,
		double	toler
	)

rsp(): diagonalize a symmetric matrix in packed (lower triangular) form in 'array'. For square symmetric matrices, see sq_rsp().

Parameters:

	nm	= rows of matrix
	n	= columns of matrix
	nv	= number of elements in lower triangle (n*(n+1)/2)
	array	= matrix to diagonalize (packed as linear array)
	e_vals	= array to hold eigenvalues
	matz	= 0 (no eigenvectors, eigenvals in ascending order) = 1 (eigenvectors and eigenvalues in ascending order) = 2 (no eigenvectors, eigenvalues in descending order) = 3 (eigenvectors and eigenvalues in descending order)
	e_vecs	= matrix of eigenvectors (one column for each eigvector)
	toler	= tolerance for eigenvalues? Often 1.0E-14.

Returns: none

Definition at line 40 of file rsp.cc.

References eigsort(), init_array(), tqli(), and tred2().

{
      int i, j, ii, ij, ierr;
      int ascend_order;
      double *fv1=NULL;
      /*double **temp=NULL;*/
      double zero = 0.0;
      double one = 1.0;
      double sw;

      /* Modified by Ed - matz can have values 0 through 3 */

      if ((matz > 3) || (matz < 0)) {
        matz = 0;
        ascend_order = 1;
        }
      else
        if (matz < 2)
          ascend_order = 1;     /* Eigenvalues in ascending order */
        else {
          matz -= 2;
          ascend_order = 0;     /* Eigenvalues in descending order */
          }

      fv1 = (double *) init_array(n);
      /*temp = (double **) init_matrix(n,n);*/

      if (n > nm) {
         ierr = 10*n;
         fprintf(stderr,"n = %d is greater than nm = %d in rsp\n",n,nm);
         exit(PSI_RETURN_FAILURE);
         }

      if (nv < (n*(n+1)/2)) {
         int num = n*(n+1)/2;
         ierr = 20*n;
         fprintf(stderr,"nv = %d is less than n*(n+1)/2 = %d in rsp\n",nv,num);
         exit(PSI_RETURN_FAILURE);
         }

      for (i=0,ij=0; i < n; i++) {
         for (j=0; j <= i; j++,ij++) {
            e_vecs[i][j] = array[ij];
            e_vecs[j][i] = array[ij];
            }
          }

      tred2(n,e_vecs,e_vals,fv1,matz);

      for (i=0; i < n; i++)
         for (j=0; j < i; j++){
            sw = e_vecs[i][j];
            e_vecs[i][j] = e_vecs[j][i];
            e_vecs[j][i] = sw;
            /*temp[i][j]=e_vecs[j][i];*/
            }
            
      tqli(n,e_vals,e_vecs,fv1,matz,toler);
      /*tqli(n,e_vals,temp,fv1,matz,toler);*/

      for (i=0; i < n; i++)
         for (j=0; j < i; j++){
            sw = e_vecs[i][j];
            e_vecs[i][j] = e_vecs[j][i];
            e_vecs[j][i] = sw;
            /*e_vecs[i][j]=temp[j][i];*/
            }

      if (ascend_order)
        eigsort(e_vals,e_vecs,n);
      else
        eigsort(e_vals,e_vecs,-n);

      free(fv1);
      /*free_matrix(temp,n);*/
}

void sq_rsp	(	int	nm,
		int	n,
		double **	array,
		double *	e_vals,
		int	matz,
		double **	e_vecs,
		double	toler
	)

sq_rsp(): diagomalize a symmetric square matrix ('array').

Parameters:

	nm	= rows of matrix
	n	= columns of matrix
	nv	= number of elements in lower triangle (n*(n+1)/2)
	array	= matrix to diagonalize
	e_vals	= array to hold eigenvalues
	matz	= 0 (no eigenvectors, eigenvals in ascending order) = 1 (eigenvectors and eigenvalues in ascending order) = 2 (no eigenvectors, eigenvalues in descending order) = 3 (eigenvectors and eigenvalues in descending order)
	e_vecs	= matrix of eigenvectors (one column for each eigvector)
	toler	= tolerance for eigenvalues? Often 1.0E-14.

Definition at line 37 of file sq_rsp.cc.

References eigsort(), free_matrix(), init_array(), init_matrix(), tqli(), and tred2().

Referenced by ael(), david(), psi::extrema::deloc::deloc(), and psi::extrema::coord_base::H_test().

{
  int i, j, ii, ij, ierr;
  int ascend_order;
  double *fv1, **temp;
  double zero = 0.0;
  double one = 1.0;

  /* Modified by Ed - matz can have the values 0 through 3 */
      
  if ((matz > 3) || (matz < 0)) {
    matz = 0;
    ascend_order = 1;
  }
  else
    if (matz < 2)
      ascend_order = 1; /* Eigenvalues in ascending order */
  else {
    matz -= 2;
    ascend_order = 0;   /* Eigenvalues in descending order */
  }

  fv1 = (double *) init_array(n);
  temp = (double **) init_matrix(n,n);

  if (n > nm) {
    ierr = 10*n;
    fprintf(stderr,"n = %d is greater than nm = %d in rsp\n",n,nm);
    exit(PSI_RETURN_FAILURE);
  }

  for (i=0; i < n; i++) {
    for (j=0; j < n; j++) {
      e_vecs[i][j] = array[i][j];
    }
  }

  tred2(n,e_vecs,e_vals,fv1,matz);

  for (i=0; i < n; i++)
    for (j=0; j < n; j++)
      temp[i][j]=e_vecs[j][i];

  tqli(n,e_vals,temp,fv1,matz,toler);

  for (i=0; i < n; i++)
    for (j=0; j < n; j++)
      e_vecs[i][j]=temp[j][i];

  if (ascend_order)
    eigsort(e_vals,e_vecs,n);
  else
    eigsort(e_vals,e_vecs,(-1)*n);

  free(fv1);
  free_matrix(temp,n);

}

void sq_to_tri	(	double **	bmat,
		double *	amat,
		int	size
	)

sq_to_tri(): converts square matrix to lower triangle

Parameters:

	bmat	= matrix to convert
	amat	= array to put lower triangle of bmat into
	size	= number of rows/columns of bmat

Returns: none

Definition at line 20 of file sq_to_tri.cc.

{
  int i, j, ij;

  ij=0;
  for(i=0; i < size; i++) {
    for(j=0 ; j <= i; j++) {
      amat[ij++] = bmat[i][j];
    }
  }

}

void tqli	(	int	n,
		double *	d,
		double **	z,
		double *	e,
		int	matz,
		double	toler
	)

tqli(): diagonalizes tridiagonal matrix output by tred2. Gives only eigenvalues if matz=0, both eigenvalues and eigenvectors if matz=1

Definition at line 21 of file tqli.cc.

Referenced by rsp(), and sq_rsp().

   {
      register int k;
      int i,j,l,m,iter;
      double dd,g,r,s,c,p,f,b,h;
      double azi;

      f=0.0;
      if (n == 1) {
         d[0]=z[0][0];
         z[0][0] = 1.0;
         return;
         }

      for (i=1; i < n ; i++) {
         e[i-1] = e[i];
         }
      e[n-1] = 0.0;
      for (l=0; l < n; l++) {
         iter = 0;
L1:
         for (m=l; m < n-1;m++) {
            dd = fabs(d[m]) + fabs(d[m+1]);
#if 0
            if (fabs(e[m])+dd == dd) goto L2;
#else
            if (fabs(e[m]) < toler) goto L2;
#endif
            }
         m=n-1;
L2:
         if (m != l) {
            if (iter++ == 30) {
               fprintf (stderr,"tqli not converging\n");
                continue;
#if 0
               exit(PSI_RETURN_FAILURE);
#endif
               }

            g = (d[l+1]-d[l])/(2.0*e[l]);
            r = sqrt(g*g + 1.0);
            g = d[m] - d[l] + e[l]/((g + DSIGN(r,g)));
            s=1.0;
            c=1.0;
            p=0.0;
            for (i=m-1; i >= l; i--) {
               f = s*e[i];
               b = c*e[i];
               if (fabs(f) >= fabs(g)) {
                  c = g/f;
                  r = sqrt(c*c + 1.0);
                  e[i+1] = f*r;
                  s=1.0/r;
                  c *= s;
                  }
               else {
                  s = f/g;
                  r = sqrt(s*s + 1.0);
                  e[i+1] = g*r;
                  c = 1.0/r;
                  s *= c;
                  }
               g = d[i+1] - p;
               r = (d[i]-g)*s + 2.0*c*b;
               p = s*r;
               d[i+1] = g+p;
               g = c*r-b;

               if (matz) {
                  double *zi = z[i];
                  double *zi1 = z[i+1];
                  for (k=n; k ; k--,zi++,zi1++) {
                     azi = *zi;
                     f = *zi1;
                     *zi1 = azi*s + c*f;
                     *zi = azi*c - s*f;
                     }
                  }
               }

            d[l] -= p;
            e[l] = g;
            e[m] = 0.0;
            goto L1;
            }
         }
   }

void tred2	(	int	n,
		double **	a,
		double *	d,
		double *	e,
		int	matz
	)

tred2(): converts symmetric matrix to a tridagonal form for use in tqli

if matz = 0, only find eigenvalues, else find both eigenvalues and eigenvectors

Returns: none

Definition at line 22 of file tred2.cc.

Referenced by rsp(), and sq_rsp().

                                                              {
    int i, j, k, l, il, ik, jk, kj;
    double f, g, h, hh, scale, scale_inv, h_inv;
    double temp;
    
    if (n == 1)
      return;
    
    for (i=n-1; i > 0; i--) {
      l = i-1;
      h = 0.0;
      scale = 0.0;
      if (l) {
        for (k=0; k <= l; k++) {
          scale += fabs(a[i][k]);
        }
        if (scale == 0.0) {
          e[i] = a[i][l];
        } else {
          scale_inv=1.0/scale;
          for (k=0; k <= l; k++) {
            a[i][k] *= scale_inv;
            h += a[i][k]*a[i][k];
          }
          f=a[i][l];
          g= -(DSIGN(sqrt(h),f));
          e[i] = scale*g;
          h -= f*g;
          a[i][l] = f-g;
          f = 0.0;
          h_inv=1.0/h;
          for (j=0; j <= l; j++) {
            if (matz)
              a[j][i] = a[i][j]*h_inv;
            g = 0.0;
            for (k=0; k <= j; k++) {
              g += a[j][k]*a[i][k];
            }
            if (l > j) {
              for (k=j+1; k <= l; k++) {
                g += a[k][j]*a[i][k];
              }
            }
            e[j] = g*h_inv;
            f += e[j]*a[i][j];
          }
          hh = f/(h+h);
          for (j=0; j <= l; j++) {
            f = a[i][j];
            g = e[j] - hh*f;
            e[j] = g;
            for (k=0; k <= j; k++) {
              a[j][k] -= (f*e[k]+ g*a[i][k]);
            }
          }
        }
      } else {
        e[i] = a[i][l];
      }
      d[i] = h;
    }
    if (matz)
      d[0] = 0.0;
    e[0] = 0.0;
    
    for (i=0; i < n; i++) {
      l = i-1;
      if (matz) {
        if (d[i]) {
          for (j=0; j <= l; j++) {
            g = 0.0;
            for (k=0; k <= l; k++) {
              g += a[i][k]*a[k][j];
            }
            for (k=0; k <= l; k++) {
              a[k][j] -= g*a[k][i];
            }
          }
        }
      }
      d[i] = a[i][i];
      if (matz) {
        a[i][i] = 1.0;
        if (l >= 0) {
          for (j=0; j<= l; j++) {
            a[i][j] = 0.0;
            a[j][i] = 0.0;
          }
        }
      }
    }
}

void tri_to_sq	(	double *	amat,
		double **	bmat,
		int	size
	)

tri_to_sq(): converts lower triangle to square matrix

Parameters:

	amat	= lower triangle matrix
	bmat	= square matrix
	size	= number of rows/cols of matrix

Returns: none

Definition at line 20 of file tri_to_sq.cc.

Referenced by psi::cscf::check_rot().

{
  int i, j, ij;

  ij=0;
  for(i = 0 ; i < size ; i++) {
    for(j = 0 ; j <= i ; j++) {
      bmat[i][j] = amat[ij];
      bmat[j][i] = amat[ij++];
    }
  }
}

void tstop

(

FILE *

outfile

)

tstop(): Stop timer

Parameters:

outfile

= output file pointer.

Definition at line 55 of file tstart.cc.

{
  int i;
  int error;
  time_t total_time;
  struct tms total_tmstime;
  char *name;
  double user_s, sys_s;

  name = (char *) malloc(40 * sizeof(char));
  error = gethostname(name, 40);
  if(error != 0) strncpy(name,"nohostname", 11);

  time_end = time(NULL);
  total_time = time_end - time_start;

  times(&total_tmstime);
  const long clk_tck = sysconf(_SC_CLK_TCK);
  user_s = ((double) total_tmstime.tms_utime)/clk_tck;
  sys_s = ((double) total_tmstime.tms_stime)/clk_tck;

  for (i=0; i < 78 ; i++) {
    fprintf(outfile,"*");
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"tstop called on %s\n", name);
  fprintf(outfile,"%s\n",ctime(&time_end));
  fprintf(outfile,"user time   = %10.2f seconds = %10.2f minutes\n",
          user_s, user_s/60.0);
  fprintf(outfile,"system time = %10.2f seconds = %10.2f minutes\n",
          sys_s, sys_s/60.0);
  fprintf(outfile,"total time  = %10d seconds = %10.2f minutes\n",
          total_time, ((double) total_time)/60.0);

  free(name);

}

void zero_arr	(	double *	a,
		int	size
	)

zero_arr(): zero out an array of length 'size'

Parameters:

	a	= array to zero out
	size	= how many elements of a to zero

Returns: none

Definition at line 21 of file zero.cc.

Referenced by pople(), and psi::detci::s1_block_fci().

{
  bzero(a,sizeof(double)*size);
}

void zero_int_array	(	int *	a,
		int	size
	)

zero_int_array() Zeroes out an array of integers 'size' integers long

Parameters:

	a	= integer array to zero out
	size	= number of elements in a to zero

Returns: none

Definition at line 59 of file int_array.cc.

Referenced by david(), ras_set(), ras_set2(), and zero_int_matrix().

{
   bzero(a,sizeof(int)*size);
}

void zero_int_matrix	(	int **	array,
		int	rows,
		int	cols
	)

zero_int_matrix(): Zero a matrix of integers. Pass the matrix, the number of rows, and the number of columns.

Parameters:

	array	= pointer to integer matrix
	rows	= number of rows in matrix
	cols	= number of columns in matrix

Returns: none

Definition at line 131 of file int_array.cc.

References zero_int_array().

{
  zero_int_array(array[0], rows*cols);
}

void zero_mat	(	double **	a,
		int	n,
		int	m
	)

zero_mat(): zero out a matrix 'a' with n rows and m columns

Parameters:

	a	= matrix of doubles to zero out
	n	= number of rows in a
	m	= number of columns in a

Definition at line 35 of file zero.cc.

Referenced by pople().

{
  register int i;

  for (i=0; i < n; i++) {
    bzero(a[i],sizeof(double)*m);
  }
}

---