---

matvec Namespace Reference

---

Detailed Description

Namespace matvec.

Compounds
class	AlleleFounderSet
	This is the "Allele Founder set" class. More...
class	AlleleOriginNode
	This is the allele origin node class. More...
class	AllelePenetranceSet
	This is the "Allele penetrance set" class. More...
class	AlleleStateNode
	This is the allele state node class. More...
class	Anterior_c
	something here More...
class	BetaDist
struct	BGNodeStruct
	This is used in SparseMatrix. More...
class	BinomialDist
class	ChiSquareDist
class	ChromStruct
	chromosome structure More...
struct	compareGNodes
	This is the comparison criterion to rank Gnodes in GNodeList based on the size of the alleleStateVector or genotypeVector. More...
struct	compareGNodesPeelId
struct	compareGNodesWeight
class	CutSet
	This is the "cut set" class. More...
class	BG
class	BGMatrix
class	Chromosome
	A chromosome. More...
class	Data
	A set of rows of data records. More...
class	DataNode
	a node in a data structure More...
class	Dblock
	This is a array of Matrix. More...
class	DisAllelePenetranceSet
	This is the "Disease Allele penetrance set" class. More...
class	DiscreteUniformDist
class	DisGenoPenetranceSet
	This is the "Disease Genotype penetrance set" class. More...
class	doubleMatrix
class	exception
class	ExponentialDist
class	FDist
class	Field
	a column for a data set More...
class	FieldStruct
	a structure for a data column More...
class	Fpmm
class	GammaDist
class	GeneticDist
	a base for Genetic distributions More...
class	GenoFounderSet
	This is the "Genotype Founder set" class. More...
class	Genome
	a genome consists of chromosomes More...
class	GenoPenetranceSet
	This is the "Genotype penetrance set" class. More...
class	GenotypeNode
	This is the genotype node class. More...
class	GeometricDist
class	GLMM
class	GNode
	This is the base graph node class. More...
struct	GNode_info
	Stores the id and distance from the cut of the GNode that will be sampled first. More...
class	GNodeList
	This is the "list of graph nodes" class. More...
class	GNodeSet
	This is the base "set of graph nodes" class. More...
class	HashNode
	a node for the hashtable More...
class	HashTable
	a hash table More...
class	IndexVector
class	Individual
	an individual consists of genomes More...
class	InvalidSample
struct	klotz0_1_
struct	klotz1_1_
class	Locus
	locus More...
exception	std
class	LocusStruct
	a structure for a locus More...
class	LogNormalDist
class	MaternalPaternalAlleles
	This is the building block class for the GenotypeNode class. More...
class	MaternalPaternalRQTLAlleles
	This is the building block class for dealing with RQTL. More...
class	Matrix
	A vector is a on-dimensional array with double precision. More...
class	Matrixwbg
class	Minimizer
class	Mixed_post_vec
class	Model
	a generalized linear mixed model More...
class	ModelTerm
	a class to hold term definition for a model More...
class	NegativeBinomialDist
struct	NodeStruct
	This is used in SparseMatrix. More...
class	NormalDist
class	NuFamily
	a nuclear family in a pedigree More...
class	Observe
struct	Parm_Elem
	It is used in Model. More...
class	ParmMap
class	Pedigree
	a pedigree More...
struct	PedNode
class	Plotter
	an object to plot curves More...
class	PoissonDist
class	Population
	a genetic population consisting of individuals More...
struct	pos_val_node
	used to store mme positions, x-values, and trait values for one obs. More...
class	Posterior_c
	something here More...
class	RAlleleFounderSet
	This is the "RAllele Founder set" class. More...
class	RAllelePenetranceSet
	This is the "RAllele penetrance set" class. More...
class	RDisAllelePenetranceSet
	This is the "RDisAllele penetrance set" class. More...
class	RecombinationSet
	This is the "Recombination set" class. More...
class	RegExp
struct	RNormal
class	RSamplerParms
class	RTransmissionSet
	This is the "RTransmission set" class. More...
class	Session
class	SparseBGMatrix
	a SparseBGMatrix is a sparse matrix of BG objects More...
class	SparseMatrix
	a sparse matrix is a matrix with very sparse elements More...
class	StatDistBase
class	Sym2x2
class	Sym3x3
class	Sym4x4
class	tDist
class	TermList
	a class to hold a list of terms in a model More...
class	TransitionSet
	This is the "Transition set" class. More...
class	TransmissionSet
	This is the "Transmission set" class. More...
class	UniformDist
class	UnknownDist
	a dummy genetic distribution More...
struct	UNormal
class	Vector
Typedefs
typedef void(*	displaytype )(const void *)
typedef double(*	penetranceType )(const Individual *, const double **)
typedef int(*	cmp_func )(const void *, const void *)
typedef int	idxtype
Functions
BG	operator+ (const BG &a, const BG &b)
BG	operator- (const BG &a, const BG &b)
BG	operator * (const BG &a, const BG &b)
BG	operator/ (double a, const BG &b)
BG	operator- (const BG &a)
BG	operator/ (const BG &a, const BG &b)
BG	sqrt (const BG &a)
BG	log (const BG &a)
BG	exp (const BG &a)
BG	operator+ (double a, const BG &b)
BG	operator- (double a, const BG &b)
BG	operator * (double a, const BG &b)
BG	operator+ (const BG &a, double b)
BG	operator- (const BG &a, double b)
BG	operator * (const BG &a, double b)
BG	operator/ (const BG &a, double b)
BG &	operator+= (BG &a, const BG &b)
BG &	operator-= (BG &a, const BG &b)
BG &	operator *= (BG &a, const BG &b)
BG &	operator/= (BG &a, const BG &b)
BG &	operator+= (BG &a, double b)
BG &	operator-= (BG &a, double b)
BG &	operator *= (BG &a, double b)
BG &	operator/= (BG &a, double b)
void	initial_BGarray (BG *v, unsigned n, double a)
void	ludcmp (BG **a, int n, int *indx, int &d, double tol)
void	lubksb (BG **a, int n, int *indx, BG *b)
bool	psdefinite (const BG **mat, const unsigned m, const unsigned n, const double tol)
Chromosome *	new_Chromosome_vec (const unsigned n)
std::ostream &	operator<< (std::ostream &stream, Data &A)
std::ostream &	operator<< (std::ostream &stream, const DataNode &A)
std::ostream &	operator<< (std::ostream &stream, Dblock &D)
int	nonsymm_eigen (const int job, double **A, const int n, double *rr, double *ri, double **vr, double **vi, const int maxiter=50)
int	symm_eigen (double *vals, double **A, const int n)
void	svdcmp (double *a[], int m, int n, double w[], double *v[])
void	ludcmp (double **a, int n, int *indx, int &d, double tol)
void	lubksb (double **a, int n, int *indx, double *b)
void	sweep (const int m, const int n, double **a, const int i0, const int i1, const double tol)
bool	psdefinite (const double **mat, const unsigned m, const unsigned n, const double tol)
Vector< bool >	operator== (const double a, const Field &b)
Vector< bool >	operator!= (const Field &a, const double b)
Vector< bool >	operator!= (const double a, const Field &b)
Vector< bool >	operator> (const Field &a, const double b)
Vector< bool >	operator> (const double a, const Field &b)
Vector< bool >	operator<= (const Field &a, const double b)
Vector< bool >	operator<= (const double a, const Field &b)
Vector< bool >	operator>= (const Field &a, const double b)
Vector< bool >	operator>= (const double a, const Field &b)
Vector< bool >	operator< (const Field &A, const double x)
Vector< bool >	operator< (const double x, const Field &A)
Field	operator+ (const Field &A, const Field &B)
Field	operator- (const Field &A, const Field &B)
Field	operator * (const Field &A, const Field &B)
Field	operator/ (const Field &A, const Field &B)
Field	operator+ (const Field &A, const double s)
Field	operator- (const Field &A, const double s)
Field	operator * (const Field &A, const double s)
Field	operator/ (const Field &A, const double s)
Field	operator+ (const double s, const Field &B)
Field	operator- (const double s, const Field &B)
Field	operator * (const double s, const Field &B)
Field	operator/ (const double s, const Field &B)
Field	sin (Field &a)
Field	asin (Field &a)
Field	cos (Field &a)
Field	acos (Field &a)
Field	tan (Field &a)
Field	atan (Field &a)
Field	ceil (Field &a)
Field	floor (Field &a)
Field	log (Field &a)
Field	log10 (Field &a)
Field	exp (Field &a)
Field	sqrt (Field &a)
Field	abs (Field &a)
Field	erf (Field &a)
Field	erfc (Field &a)
int	compare_DataNode (const void *A, const void *B)
std::ostream &	operator<< (std::ostream &stream, const Field &V)
int	fieldindex (const std::string &fw, const FieldStruct *field_vec, const unsigned n)
Genome *	new_Genome_vec (const unsigned m, GeneticDist *D)
unsigned	ginverse1 (double **a, const unsigned n, double &lgdet, int mode, const double tol)
unsigned	ginverse1 (BG **a, const unsigned n, BG &lgdet, int mode, const double tol)
void	sweep (const int m, const int n, BG **a, const int i0, const int i1, const BG tol)
doubleMatrix *	KP (Data *D, doubleMatrix &SKM, const std::string &lpl, const std::string &censor)
doubleMatrix	operator * (Vector< double > &a, Vector< double > &b)
void	bound_pd (double *varnew, int ntrait, doubleMatrix &P)
void	getlambda (double **lambda, const int n)
void	normal_loginClasses (matvec::Observe &observe)
int	squeeze (unsigned len, int *ia, int *ja, double *a)
void	hsort_ija (unsigned n, int *ia, int *ja, double *a)
long	next_prime (long n)
HashNode *	new_HashNode_vec (const unsigned n, const size_t s)
double	ranf (void)
int	compare_ind_pt (const void *x, const void *y)
int	compare_ind_id (const void *a, const void *b)
int	compare_ind_gid (const void *a, const void *b)
int	compare_mother_id (const void *x, const void *y)
int	compare_father_id (const void *x, const void *y)
RNormal *	new_n_posterior_vec (const int m)
template<class T> Matrix< T >	operator/ (const T &a, const Matrix< T > &b)
template<class T> Vector< T >	operator * (const Vector< T > &a, const Matrix< T > &b)
template<class T> Matrix< T >	hadjoin (const Matrix< T > &a, const Matrix< T > &b)
template<class T> Matrix< T >	hadjoin (const Matrix< T > &a, const Vector< T > &b)
template<class T> Matrix< T >	hadjoin (const Vector< T > &a, const Matrix< T > &b)
template<class T> Matrix< T >	vadjoin (const Matrix< T > &a, const Matrix< T > &b)
template<class T> Matrix< T >	vadjoin (const Matrix< T > &a, const Vector< T > &b)
template<class T> Matrix< T >	vadjoin (const Vector< T > &a, const Matrix< T > &b)
template<class T> Matrix< T >	diag (const Vector< T > &a)
template<class T> Matrix< T >	diag (const Matrix< T > &a)
template<class T> Matrix< T >	sin (const Matrix< T > &a)
template<class T> Matrix< T >	asin (const Matrix< T > &a)
template<class T> Matrix< T >	cos (const Matrix< T > &a)
template<class T> Matrix< T >	acos (const Matrix< T > &a)
template<class T> Matrix< T >	tan (const Matrix< T > &a)
template<class T> Matrix< T >	atan (const Matrix< T > &a)
template<class T> Matrix< T >	ceil (const Matrix< T > &a)
template<class T> Matrix< T >	floor (const Matrix< T > &a)
template<class T> Matrix< T >	log (const Matrix< T > &a)
template<class T> Matrix< T >	log10 (const Matrix< T > &a)
template<class T> Matrix< T >	exp (const Matrix< T > &a)
template<class T> Matrix< T >	sqrt (const Matrix< T > &a)
template<class T> Matrix< T >	abs (const Matrix< T > &a)
template<class T> Matrix< T >	erf (const Matrix< T > &a)
template<class T> Matrix< T >	erfc (const Matrix< T > &a)
template<class T> bool	all (const Matrix< T > &a)
template<class T> bool	any (const Matrix< T > &a)
template<class T> Matrix< T >	operator+ (const T &a, const Matrix< T > &b)
template<class T> Matrix< T >	operator- (const T &a, const Matrix< T > &b)
template<class T> Matrix< T >	operator * (const T &a, const Matrix< T > &b)
template<class T> std::ostream &	operator<< (std::ostream &os, const Matrix< T > &a)
void	linmin (double *p, double *xi, int n, double &fret)
double	f1dim (double x)
double	brent (double ax, double bx, double cx, double(*f)(double), double tol, double &xmin)
void	mnbrak (double &ax, double &bx, double &cx, double &fa, double &fb, double &fc, double(*func)(double))
double	funk (Vector< double > &x)
double	nr_powell (Data *D, Model *M, Vector< double > &p, Matrix< double > &xi, int n, int &iter, double ftol)
void	svd (const int n, double eps, const double tol, double **ab, double *q)
void	sort (const int n, double *d, double **v)
void	print (const int n, const Vector< double > &x)
void	flin (const double l, const int j, const Vector< double > &x, const int n, const Matrix< double > &v, Vector< double > &w, const Matrix< double > &q01)
void	nest_xaction (std::string &s)
unsigned	lidx (const unsigned i, const unsigned j, const unsigned n)
void	ginverse2 (double *a, const unsigned n, const double tol)
void	preconditioning (double *result, const double **M, const double *r, const unsigned n)
void	ABCt (double **A, const unsigned ma, const unsigned na, double **B, const unsigned mb, const unsigned nb, double **C, const unsigned mc, const unsigned nc, double **OUT, const unsigned mo, const unsigned no)
void	balance (double **mat, const int n, int *low, int *hi, double *scale, const double base)
void	elemhess (int job, double **mat, const int n, const int low, const int hi, double **vr, double **vi)
void	unbalance (const int n, double **vr, double **vi, const int low, const int hi, double *scale)
int	nonsymm_eigen_core (const int job, double **mat, const int n, const int low, const int hi, double *valr, double *vali, double **vr, double **vi, const double eps, const int maxiter)
int	compare_ind_ptr (const void *a, const void *b)
NuFamily *	new_NuFamily_vec (const unsigned m)
int	pedcp1 (const void *a, const void *b)
int	pedcp2 (const void *a, const void *b)
std::ostream &	operator<< (std::ostream &stream, Pedigree &A)
unsigned	sampling (double *dist_prob, const unsigned size)
void	G_to_gamete (const unsigned G_id, unsigned k1, unsigned k2, unsigned size)
int	pdm_val (const unsigned *goff, const unsigned *gsire, const unsigned *gdam, Matrix< double > &pdm, const int ori)
void	compute_pdq (const double er, const Matrix< double > &pdm, Matrix< double > &pdq)
double	pr_osd (const unsigned goff[], const unsigned gsire[], const unsigned gdam[], const int ori)
int	inva22 (Matrix< double > &a22)
double	penetrance (const Individual *I, const double **res_var)
int	compare_seq_id (const void *a, const void *b)
int	compare_nuf_ptr (const void *a, const void *b)
int	comp1 (const int *i, const int *j)
void	gs_ioc (int *ia, int *ja, BG *a, const BG *rhs, BG *sol, const int n, const double relax, const double stopval, const int mxiter, const double tol, unsigned &cov)
int	squeeze (unsigned len, int *iiv, int *jjv, BG *aav)
void	hsort_ija (unsigned n, int *iiv, int *jjv, BG *aav)
int	operator== (const SparseBGMatrix &A, const SparseBGMatrix &B)
int	operator!= (const SparseBGMatrix &A, const SparseBGMatrix &B)
Vector< BG >	operator * (const SparseBGMatrix &A, const Vector< BG > &v)
SparseBGMatrix	operator * (const SparseBGMatrix &A, const BG s)
SparseBGMatrix	operator * (const BG s, const SparseBGMatrix &A)
void	sparse_dense (unsigned n, int *ia, int *ja, BG *a, BG **matx)
int	comp (const int *i, const int *j)
void	gs_ioc (int *ia, int *ja, double *a, const double *rhs, double *sol, const int n, const double relax, const double stopval, const int mxiter, const double tol, unsigned &cov)
int	operator== (const SparseMatrix &A, const SparseMatrix &B)
int	operator!= (const SparseMatrix &A, const SparseMatrix &B)
Vector< double >	operator * (const SparseMatrix &A, const Vector< double > &v)
SparseMatrix	operator * (const SparseMatrix &A, const double s)
SparseMatrix	operator * (const double s, const SparseMatrix &A)
void	sparse_dense (unsigned n, int *ia, int *ja, double *a, double **matx)
double	snorm (void)
double	sexpo (void)
double	sgamma (const double a)
double	genbet (const double aa, const double bb)
double	genchi (const int df)
double	genexp (const double av)
double	genf (const int dfn, const int dfd)
double	gengam (const double a, const double r)
double	gennch (const int df, const double xnonc)
double	gennf (const int dfn, const int dfd, const double xnonc)
double	gennor (const double av, const double sd)
double	genunf (const double low, const double high)
long	ignbin (const long n, const double pp)
long	ignpoi (const double mu)
long	ignuin (const long low, const long high)
void	set_seed (const unsigned iseed)
double	ran1 (void)
double	gasdev (void)
double	gamdev (const double A, const double R)
double	factrl (const long n)
double	BinCoef (const long n, const long k)
double	gammln (const double xx)
double	gammin (const double x, const double a)
double	betain (const double x, const double a, const double b, const double beta)
double	Normal_cdf (const double x)
double	Gamma_cdf (const double g, const double alfa, const double theta)
double	Beta_cdf (const double x, const double alfa, const double beta, const double lambda, int &errcode)
double	F_cdf (const double f, const double df1, const double df2, const double nc, int &errcode)
double	ChiSquare_cdf (const double x, const double df, const double theta, int &errcode)
double	t_cdf (const double t, const double df, const double delta)
double	t_expected_value (const double df, const double nc)
double	Normal_inv (const double p)
double	ChiSquare_inv (const double p, const double df, const double nc)
double	t_inv (const double p, const double df, const double nc)
double	F_inv (const double p, const double df1, const double df2, const double nc)
void	warning (const char format[],...)
double	betacf (const double a, const double b, const double x)
int	zufall_ (int n, double *a)
int	zufalli_ (int seed)
double	fsign (double num, double sign)
void	tred2 (double **a, int n, double *d, double *e)
void	tqli (double *d, double *e, int n, double **z)
void	check_ptr (void *a)
int	split (const std::string &str, const std::string &sep, std::vector< std::string > *strvec)
bool	validline (char *line)
unsigned	est_nze (const unsigned dim, const unsigned ntrait, const int pedsize)
Variables
int	DISCRETE_TRAIT
doubleMatrix	PENETRANCE_MATRIX
double	at
double	bt
double	ct
double	maxarg1
double	maxarg2
doubleMatrix *	P_Mat = NULL
double	global_mu
double	global_vg
double	global_ve
Data *	myD
Model *	myM
int	nfunk
double	sqrarg
int	ncom = 0
double *	pcom = 0
double *	xicom = 0
double	qd [2]
double	fx
double	qf1
double	macheps
double	t
double	m2
double	m4
double	SmaLL
double	vsmall
double	vlarge
double	large
double	dmin
double	ldt
double	h
int	nf
int	nl
double	MIXED_TOL = 0.00000000001
int	INPUT_LINE_WIDTH = 1024
double	MY_NEG_HUGE = -1.0e+20
Session	SESSION
struct {
int   fill_1 [1214]
int   e_2
}	klotz0_
struct {
double   fill_1 [1024]
int   e_2 [2]
double   e_3
}	klotz1_
long	RAND_STATE_ARRAY [64]
const int	MATVEC_MESSAGE_BUF_LEN = 1024

---

Typedef Documentation

typedef int(* matvec::cmp_func)(const void *, const void *)

Definition at line 43 of file sparsebgmatrix.h. Referenced by matvec::SparseMatrix::minfilsymelm_node_list(), and matvec::SparseBGMatrix::minfilsymelm_node_list().

typedef void(* matvec::displaytype)(const void*)

Definition at line 60 of file hashtable.h. Referenced by matvec::HashTable::display().

typedef int matvec::idxtype

Definition at line 55 of file sparsebgmatrix.h. Referenced by matvec::SparseMatrix::minfilsymelm_node_list(), and matvec::SparseBGMatrix::minfilsymelm_node_list().

typedef double(* matvec::penetranceType)(const Individual*, const double **)

Definition at line 40 of file individual.h. Referenced by matvec::Population::setPenetrance(), and matvec::Individual::setPenetrance().

---

Function Documentation

void ABCt (  double **    A, const unsigned    ma, const unsigned    na, double **    B, const unsigned    mb, const unsigned    nb, double **    C, const unsigned    mc, const unsigned    nc, double **    OUT, const unsigned    mo, const unsigned    no )

template<class T> Matrix<T> abs (  const Matrix< T > &    a )

Definition at line 1114 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::abs); }

Field abs (  Field &    a )

Definition at line 774 of file field.cpp. References matvec::Field::map(). Referenced by matvec::doubleMatrix::norm(). {return a.map(std::fabs);}

template<class T> Matrix<T> acos (  const Matrix< T > &    a )

Definition at line 1105 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::acos); }

Field acos (  Field &    a )

Definition at line 765 of file field.cpp. References matvec::Field::map(). {return a.map(std::acos);}

template<class T> bool all (  const Matrix< T > &    a )

Definition at line 1118 of file matrix.h. References matvec::Matrix< T >::all(). { return a.all();}

template<class T> bool any (  const Matrix< T > &    a )

Definition at line 1119 of file matrix.h. References matvec::Matrix< T >::any(). { return a.any();}

template<class T> Matrix<T> asin (  const Matrix< T > &    a )

Definition at line 1103 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::asin); }

Field asin (  Field &    a )

Definition at line 763 of file field.cpp. References matvec::Field::map(). {return a.map(std::asin);}

template<class T> Matrix<T> atan (  const Matrix< T > &    a )

Definition at line 1107 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::atan); }

Field atan (  Field &    a )

Definition at line 767 of file field.cpp. References matvec::Field::map(). {return a.map(std::atan);}

void balance (  double **    mat, const int    n, int *    low, int *    hi, double *    scale, const double    base )

Definition at line 28 of file nonsymm_eigen.cpp. Referenced by nonsymm_eigen(). { //////////////////////////////////////////////////////////////////// // Balance a matrix for calculation of eigenvalues and eigenvectors // //////////////////////////////////////////////////////////////////// double c,f,g,r,s; int i,j,k,l,done; /* search for rows isolating an eigenvalue and push them down */ for (k = n - 1; k >= 0; k--) { for (j = k; j >= 0; j--) { for (i = 0; i <= k; i++) { if (i != j && fabs(mat[j][i]) != 0) break; } if (i > k) { scale[k] = j; if (j != k) { for (i = 0; i <= k; i++) { c = mat[i][j]; mat[i][j] = mat[i][k]; mat[i][k] = c; } for (i = 0; i < n; i++) { c = mat[j][i]; mat[j][i] = mat[k][i]; mat[k][i] = c; } } break; } } if (j < 0) break; } /* search for columns isolating an eigenvalue and push them left */ for (l = 0; l <= k; l++) { for (j = l; j <= k; j++) { for (i = l; i <= k; i++) { if (i != j && fabs(mat[i][j]) != 0) break; } if (i > k) { scale[l] = j; if (j != l) { for (i = 0; i <= k; i++) { c = mat[i][j]; mat[i][j] = mat[i][l]; mat[i][l] = c; } for (i = l; i < n; i++) { c = mat[j][i]; mat[j][i] = mat[l][i]; mat[l][i] = c; } } break; } } if (j > k) break; } *hi = k; *low = l; /* balance the submatrix in rows l through k */ for (i = l; i <= k; i++) scale[i] = 1; do { for (done = 1,i = l; i <= k; i++) { for (c = 0,r = 0,j = l; j <= k; j++) { if (j != i) { c += fabs(mat[j][i]); r += fabs(mat[i][j]); } } if (c != 0 && r != 0) { g = r / base; f = 1; s = c + r; while (c < g) { f *= base; c *= base * base; } g = r * base; while (c >= g) { f /= base; c /= base * base; } if ((c + r) / f < 0.95 * s) { done = 0; g = 1 / f; scale[i] *= f; for (j = l; j < n; j++) mat[i][j] *= g; for (j = 0; j <= k; j++) mat[j][i] *= f; } } } } while (!done); }

double matvec::Beta_cdf (  const double    x, const double    alfa, const double    beta, const double    lambda, int &    errcode )

Definition at line 491 of file stat1.cpp. References betain(), ERRMAX, exp(), gammln(), ITMAX, log(), sqrt(), and warning(). Referenced by matvec::BetaDist::cdf(), and F_cdf(). { //////////////////////////////////////////////////////////////////// // ALGORITHM AS226 APPL. STATIST. (1987) VOL. 36, NO. 2 // Incorporates modification AS R84 from AS vol. 39, pp311-2, 1990 // // Returns the cumulative probability of X for the non-central beta // distribution with parameters A, B and non-centrality LAMBDA // // Auxiliary routines required: ALOGAM - log-gamma function (ACM // 291 or AS 245), and BETAIN - incomplete-beta function (AS 63) // // errcode = 0 no error // = 1 parameter range error // = 2 fail to converge ///////////////////////////////////////////////////////////////////// if (alfa <= 0.0 || beta <= 0.0) throw exception(" Beta_cdf(): bad args, value > 0 expected"); errcode = 0; if (x <= 0.0) return 0.0; if (x >= 1.0) return 1.0; double lnbeta, betanc = 0.0; if (lambda == 0.0) { lnbeta = gammln(alfa) + gammln(beta) - gammln(alfa+beta); betanc = betain(x,alfa,beta,lnbeta); return betanc; } double c,x0,gx,q,xj,ax,sumq,temp,alfa0,errbd; c = lambda*0.5; x0 = 0.0; // Initialize the series ... q = c - 5.0*sqrt(c); if (q > 0.0) x0 = static_cast<int>(q); alfa0 = alfa + x0; lnbeta = gammln(alfa0) + gammln(beta) - gammln(alfa0 + beta); temp = betain(x,alfa0,beta,lnbeta); gx = exp(alfa0*log(x) + beta*log(1.0 - x) - lnbeta - log(alfa0)); if (alfa0 > alfa) { q = exp(-c + x0*log(c) - gammln(x0 + 1.0)); } else { q = exp(-c); } xj = 0.0; ax = q*temp; sumq = 1.0 - q; betanc = ax; do { // recur over subsequent terms until convergence is achieved... xj++; temp -= gx; gx *= x*(alfa+ beta + xj - 1.0)/(alfa + xj); q *= c/xj; sumq -= q; ax = temp*q; betanc += ax; errbd = (temp - gx)*sumq; } while(static_cast<int>(xj) < ITMAX && errbd > ERRMAX); if (errbd > ERRMAX) { warning("Beta_cdf(): fail to converge"); errcode = 2; } if (betanc > 1.0) betanc = 1; return betanc; }

double betacf (  const double    a, const double    b, const double    x )

Definition at line 176 of file stat1.cpp. References matvec::Session::epsilon, ITMAX, SESSION, and warning(). Referenced by betain(). { double qap,qam,qab,em,tem,d; double bz,bm=1.0,bp,bpp; double az=1.0,am=1.0,ap,app,aold; int m; qab=a+b; qap=a+1.0; qam=a-1.0; bz=1.0-qab*x/qap; for (m=1; m<=ITMAX; m++) { em = static_cast<double>(m); tem = em+em; d = em*(b-em)*x/((qam+tem)*(a+tem)); ap = az+d*am; bp = bz+d*bm; d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem)); app = ap+d*az; bpp = bp+d*bz; aold = az; am = ap/bpp; bm = bp/bpp; az = app/bpp; bz = 1.0; if (fabs(az-aold) < (10.0*SESSION.epsilon*fabs(az))) return az; } warning("betacf(), not converged"); abort(); return 0.0; }

double matvec::betain (  const double    x, const double    a, const double    b, const double    beta )

Definition at line 290 of file stat1.cpp. References betacf(), matvec::Session::epsilon, exp(), log(), and SESSION. Referenced by Beta_cdf(), matvec::BinomialDist::sample(), and t_cdf(). { ////////////////////////////////////////////////////////// // return incomplete beta function for arguments: // x between 0 and 1, // a and b positive, // beta, log of complete beta function value /////////////////////////////////////////////////////////// if (x < 0.0 || x > 1.0) throw exception(" betain(%f): bad arg1, value between (0,1) expected"); if (x == 0.0) return 0.0; if (x == 1.0) return 1.0; double fv,bt; if (x < SESSION.epsilon || fabs(x - 1.0) < SESSION.epsilon) { bt=0.0; } else { fv = a*log(x) + b*log(1.0-x) - beta; if (fv < -300) { bt = 0.0;} else { bt = exp(fv);} } if (x < (a+1.0)/(a+b+1.0)) { fv = bt*betacf(a,b,x)/a; } else { fv = 1.0-bt*betacf(b,a,1.0-x)/b; } return fv; }

double matvec::BinCoef (  const long    n, const long    k )

Definition at line 328 of file stat1.cpp. References gammln(). { if (n<0L || k<0L || k > n) throw exception("BinCoef(): bad args"); double a,b,c; a = gammln(static_cast<double>(n+1)); b = gammln(static_cast<double>(k+1)); c = gammln(static_cast<double>(n-k)); return std::floor(0.5 + std::exp(a - b - c)); }

void matvec::bound_pd (  double *    varnew, int    ntrait, doubleMatrix &    P )

Definition at line 1276 of file glmm2.cpp. References diag(), matvec::doubleMatrix::eigen(), matvec::Vector< T >::resize(), and matvec::Matrix< double >::transpose(). { int i, j, k, t1, t2; doubleMatrix varbound (numtrait, numtrait); Vector<double> eig, flg; flg.resize (numtrait); for (k = 0, i = 0; i < numtrait; i++) for (j = i; j < numtrait; j++, k++) { varbound[i][j] = varnew[k]; varbound[j][i] = varnew[k]; } P = varbound; eig = P.eigen (); for (t1 = 0; t1 < numtrait; t1++) { flg[t1] = 0; if (eig[t1] < 1.e-3) { eig[t1] = .99e-3; flg[t1] = 1; } if (eig[t1] > 1.e10) { eig[t1] = 1.01e10; flg[t1] = 0; } } varbound = P * diag(eig) * P.transpose(); for (t1 = 0; t1 < numtrait; t1++) { if (flg[t1]) { for (t2 = 0; t2 < numtrait; t2++) { P[t1][t2] = 0; } } } for (k = 0, t1 = 0; t1 < numtrait; t1++) for (t2 = t1; t2 < numtrait; t2++, k++) { varnew[k] = varbound[t1][t2]; } }

double matvec::brent (  double    ax, double    bx, double    cx, double(*    f)(double), double    tol, double &    xmin )  [static]

Definition at line 248 of file min_nr_powell.cpp. References CGOLD, ITMAX, SHFT, SIGN, warning(), and ZEPS. Referenced by linmin(). { int iter; double a,b,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm; double e=0.0; double d=0.0; a=((ax < cx) ? ax : cx); b=((ax > cx) ? ax : cx); x=w=v=bx; fw=fv=fx=(*f)(x); for (iter=1;iter<=ITMAX;iter++) { xm=0.5*(a+b); tol2=2.0*(tol1=tol*fabs(x)+ZEPS); if (fabs(x-xm) <= (tol2-0.5*(b-a))) { xmin=x; return fx; } if (fabs(e) > tol1) { r=(x-w)*(fx-fv); q=(x-v)*(fx-fw); p=(x-v)*q-(x-w)*r; q=2.0*(q-r); if (q > 0.0) p = -p; q=fabs(q); etemp=e; e=d; if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) d=CGOLD*(e=(x >= xm ? a-x : b-x)); else { d=p/q; u=x+d; if (u-a < tol2 || b-u < tol2) d=SIGN(tol1,xm-x); } } else { d=CGOLD*(e=(x >= xm ? a-x : b-x)); } u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d)); fu=(*f)(u); if (fu <= fx) { if (u >= x) { a=x; } else { b=x; } SHFT(v,w,x,u) SHFT(fv,fw,fx,fu) } else { if (u < x) { a=u; } else { b=u; } if (fu <= fw || w == x) { v=w; w=u; fv=fw; fw=fu; } else if (fu <= fv || v == x || v == w) { v=u; fv=fu; } } } warning(" Too many iterations in BRENT"); xmin=x; return fx; }

template<class T> Matrix<T> ceil (  const Matrix< T > &    a )

Definition at line 1108 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::ceil); }

Field ceil (  Field &    a )

Definition at line 768 of file field.cpp. References matvec::Field::map(). {return a.map(std::ceil);}

void matvec::check_ptr (  void *    a )

Definition at line 37 of file util.cpp. Referenced by matvec::Model::build_model_term(), matvec::Population::build_nufamily(), matvec::Population::build_trans_mat(), matvec::RegExp::copy(), matvec::Model::copyfrom(), elemhess(), matvec::RegExp::find(), matvec::Vector< T >::initialize(), matvec::Pedigree::input(), matvec::Pedigree::input_pedsheet(), matvec::RegExp::match(), matvec::HashTable::maxsize(), matvec::GeneticDist::nchrom(), new_Genome_vec(), matvec::Data::newcol(), matvec::Population::peeling_sequence(), matvec::Model::prepare_data(), matvec::Model::prior_dist(), matvec::Model::re_hash_data(), matvec::RegExp::replace(), matvec::Vector< T >::reserve(), matvec::Vector< T >::resize(), matvec::TermList::resize(), matvec::HashTable::resize(), matvec::Data::resize(), matvec::Model::RSamplerPrior_dist(), matvec::Population::sub(), and matvec::Model::variance(). { ; }

double matvec::ChiSquare_cdf (  const double    x, const double    df, const double    theta, int &    errcode )

Definition at line 360 of file stat1.cpp. References ERRMAX, exp(), gammin(), gammln(), ITRMAX, and log(). Referenced by matvec::ChiSquareDist::cdf(), ChiSquare_inv(), and matvec::GLMM::contrast(). { /////////////////////////////////////////////////////////////////////// // ALGORITHM AS 275 APPL.STATIST. (1992), VOL.41, NO.2 // // Computes the noncentral chi-square cumulative distribution function // with positive real degrees of freedom df and nonnegative // noncentrality parameter theta // // errcode = 0 no error // = 1 parameter range error // = 2 fail to converge // Tianlin Wang, UIUC, Mon Jun 6 14:02:11 CDT 1994 //////////////////////////////////////////////////////////////////////// if (df <= 0.0) throw exception(" ChiSquare_cdf(): bad arg2: %f, value > 0 expected"); if (theta < 0.0) throw exception(" ChiSquare_cdf(): bad arg3, nonnegative value expected"); errcode = 0; if (x <= 0.0) return 0.0; double chi2nc = x; if (theta == 0.0) return gammin(chi2nc*0.5, df*0.5); double lam,u,v,x2,f2,t,errbd,n2; lam = 0.5*theta; unsigned n = 1; // Evaluate the first term u = exp(-lam); v = u; x2 = 0.5*x; f2 = 0.5*df; t = f2*log(x2) - x2 - gammln(f2+1.0); if (t <= -600.0) return 1.0; // this is added in case exp(t) == 0 t = exp(t); chi2nc = v*t; n2 = static_cast<double>(n + n); while ((df + n2) <= x) { u *= lam/static_cast<double>(n); // Evaluate the next term of the expansion v += u; t *= x/(df + n2); chi2nc += v*t; n++; n2 = static_cast<double>(n + n); } do { n2 = static_cast<double>(n + n); errbd = t*x/(df + n2 -x); if (n>ITRMAX) throw exception ("ChiSquare_cdf(): fail to converge"); u *= lam/static_cast<double>(n); // Evaluate the next term of the expansion v += u; t *= x/(df + n2); chi2nc += v*t; n++; } while (errbd > ERRMAX); return chi2nc; }

double matvec::ChiSquare_inv (  const double    p, const double    df, const double    nc )

Definition at line 701 of file stat1.cpp. References ChiSquare_cdf(), exp(), gammin(), gammln(), log(), Normal_inv(), and sqrt(). Referenced by matvec::ChiSquareDist::inv(). { ////////////////////////////////////////////////////////////////// // // Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35 // // To evaluate the percentage points of the chi-squared // probability distribution function. // // Incorporates the suggested changes in AS R85 (vol.40(1), // pp.233-5, 1991) // // Auxiliary routines required: PPND = AS 111 (or AS 241) and // GAMMAD = AS 239. // NOTE: // ChiSquare_inv(1.0e-16,100,0.0) returns 52.4813, it seems not accurate /////////////////////////////////////////////////////////////////////// if (p < 0.0 || p >= 1.0 ) throw exception("ChiSquare_inv(): bad arg1, value in (0,1) expected"); if (df <= 0.0) throw exception(" ChiSquare_inv(): bad arg2, value > 0 expected"); if (p == 0.0) return 0.0; int errcode = 0; double ppchi2 = 0.0; if (nc != 0.0) { double eps = 1.0e-7; double bot = 0.0; double top = 26000.0; double p2; while (top-bot > eps) { ppchi2 = (top + bot)*0.5; p2 = ChiSquare_cdf(ppchi2,df,nc,errcode); if (errcode != 0) break; if (p2 < p ) { bot = ppchi2; } else if (p2 > p) { top = ppchi2; } else { break; } } return ppchi2; } double one = 1.0; double aa,e,a,b,d,ch,p1,p2,q,s1,s2,s3,s4,s5,s6,t,x,xx; aa = 0.6931471806e+0; e = 0.5e-6; double* c = new double [38]; c[0]= 0.01; c[1]= 0.222222; c[2]= 0.32; c[3]= 0.4; c[4]= 1.24; c[5]= 2.2; c[6]= 4.67; c[7]= 6.66; c[8]= 6.73; c[9]= 13.32; c[10]= 60.0; c[11]= 70.0; c[12]= 84.0; c[13]= 105.0; c[14]= 120.0; c[15]= 127.0; c[16]= 140.0; c[17]= 1175.0; c[18]= 210.0; c[19]= 252.0; c[20]= 2264.0; c[21]= 294.0; c[22]= 346.0; c[23]= 420.0; c[24]= 462.0; c[25]= 606.0; c[26]= 672.0; c[27]= 707.0; c[28]= 735.0; c[29]= 889.0; c[30]= 932.0; c[31]= 966.0; c[32]= 1141.0; c[33]= 1182.0; c[34]= 1278.0; c[35]= 1740.0; c[36]= 2520.0; c[37]= 5040.0; xx = df*0.5; double g = gammln(xx); d = xx - one; // starting approximation for small chi-squared if (df < -c[4]*log(p)) { ch = pow(p*xx*exp(g + xx*aa),one/xx); if (ch <= e) { if(c){ delete [] c; c=0; } return ch;} } else { if (df < c[2]) { // starting approximation for df <= 0.32 ch = c[3]; a = log(one-p); do { q = ch; p1 = one + ch*(c[6] + ch); p2 = ch*(c[8] + ch*(c[7] + ch)); t = (c[6] + 2.0*ch)/p1 - (c[8] + ch*(c[9] + 3.0*ch))/p2 - 0.5; ch -= (1.0 - exp(a + g + 0.5*ch + d*aa)*p2/p1)/t; } while (fabs(q/ch - one) > c[0]); } else { // call to algorithm AS 111 - note that p has been tested above. // AS 241 could be used as an alternative. x = Normal_inv(p); // starting approximation using Wilson and Hilferty estimate p1 = c[1]/df; ch = df*pow(x*sqrt(p1) + one - p1,3.0); // starting approximation for p tending to 1 if (ch > c[5]*df + 6.0) ch = -2.0*(log(one-p) - d*log(ch*0.5) + g); } } int i = 0; while (i++ < 20) { // call to algorithm AS 239 and calculation of seven term Taylor series q = ch; p1 = 0.5*ch; p2 = p - gammin(p1,xx); t = p2*exp(xx*aa + g + p1 - d*log(ch)); b = t/ch; a = 0.5*t - b*d; s1 = (c[18]+ a*(c[16]+ a*(c[13]+ a*(c[12]+ a*(c[11]+ c[10]*a)))))/c[23]; s2 = (c[23]+ a*(c[28]+ a*(c[31]+ a*(c[32] + c[34]*a))))/c[36]; s3 = (c[18]+ a*(c[24]+ a*(c[27]+ c[30]*a)))/c[36]; s4 = (c[19]+ a*(c[26]+ c[33]*a)+ d*(c[21]+ a*(c[29]+ c[35]*a)))/c[37]; s5 = (c[12]+ c[20]*a+ d*(c[17]+ c[25]*a))/c[36]; s6 = (c[14]+ d*(c[22]+ c[15]*d))/c[37]; ch += t*(one + 0.5*t*s1- b*d*(s1- b*(s2- b*(s3- b*(s4- b*(s5- b*s6)))))); if (fabs(q/ch - one) > e) break; } ppchi2 = ch; if(c){ delete [] c; c=0; } return ppchi2; }

int matvec::comp (  const int *    i, const int *    j )

Definition at line 33 of file sparsematrix.cpp. Referenced by matvec::SparseMatrix::minfilsymelm_node_list(). { return *i-*j; }

int matvec::comp1 (  const int *    i, const int *    j )

Definition at line 33 of file sparsebgmatrix.cpp. Referenced by matvec::SparseBGMatrix::minfilsymelm_node_list(). { return *i-*j; }

int compare_DataNode (  const void *    A, const void *    B )

Definition at line 969 of file field.cpp. References matvec::DataNode::double_val(), and matvec::DataNode::missing. Referenced by matvec::Field::sort(). { DataNode *x = (DataNode *)A; DataNode *y = (DataNode *)B; if (x->missing && y->missing) { return 0; } else if (x->missing && !y->missing) { return 1; } else if (!x->missing && y->missing) { return -1; } else { if (x->double_val() < y->double_val()) { return -1; } else if (x->double_val() == y->double_val()) { return 0; } else { return 1; } } }

int matvec::compare_father_id (  const void *    x, const void *    y )

Definition at line 80 of file individual.cpp. Referenced by matvec::Population::build_spouse_info(). { Individual **U1,**U2; U1 = (Individual **)x; U2 = (Individual **)y; return ((*U1)->father_id() - (*U2)->father_id()); }

int matvec::compare_ind_gid (  const void *    a, const void *    b )

Definition at line 63 of file individual.cpp. Referenced by matvec::Population::renum(). { Individual **x,**y; x = (Individual **)a; y = (Individual **)b; return ((*x)->gid() - (*y)->gid()); }

int matvec::compare_ind_id (  const void *    a, const void *    b )

Definition at line 54 of file individual.cpp. { Individual **x,**y; x = (Individual **)a; y = (Individual **)b; return ((*x)->id() - (*y)->id()); }

int matvec::compare_ind_pt (  const void *    x, const void *    y )

Definition at line 40 of file individual.cpp. { Individual **U1,**U2; U1 = (Individual **)x; U2 = (Individual **)y; int retval = 1; if (*U1 < *U2) retval = -1; else if (*U1 == *U2) { retval = 0; } return retval; }

int compare_ind_ptr (  const void *    a, const void *    b )  [static]

Definition at line 51 of file nufamily.cpp. { Individual **x,**y; x = (Individual **)a; y = (Individual **)b; if (x[0] == y[0]) { return 0; } else { return 1; } }

int matvec::compare_mother_id (  const void *    x, const void *    y )

Definition at line 72 of file individual.cpp. Referenced by matvec::Population::build_spouse_info(). { Individual **U1,**U2; U1 = (Individual **)x; U2 = (Individual **)y; return ((*U1)->mother_id() - (*U2)->mother_id()); }

int compare_nuf_ptr (  const void *    a, const void *    b )  [static]

Definition at line 39 of file pop_peeling.cpp. Referenced by matvec::Population::break_loop(). { NuFamily **x,**y; x = (NuFamily **)a; y = (NuFamily **)b; if (x[0] == y[0]) return 0; else return 1; }

int compare_seq_id (  const void *    a, const void *    b )

Definition at line 32 of file pop_peeling.cpp. References matvec::NuFamily::seq_id. Referenced by matvec::Population::detect_loop(). { NuFamily **x,**y; x = (NuFamily **)a; y = (NuFamily **)b; return (x[0]->seq_id - y[0]->seq_id); }

void matvec::compute_pdq (  const double    er, const Matrix< double > &    pdm, Matrix< double > &    pdq )

Definition at line 831 of file pop_mblup.cpp. Referenced by matvec::Population::add_Gv_inv(), matvec::Population::relv_inv(), and matvec::Population::relv_nomissing(). { // C = S*R, where C is the PDQ matrix of 2 by 4 for (int i=0; i<2; i++) { pdq[i][0] = (1-er)*pdm[i][0] + er *pdm[i][1]; pdq[i][1] = er *pdm[i][0] + (1-er)*pdm[i][1]; pdq[i][2] = (1-er)*pdm[i][2] + er *pdm[i][3]; pdq[i][3] = er *pdm[i][2] + (1-er)*pdm[i][3]; } }

template<class T> Matrix<T> cos (  const Matrix< T > &    a )

Definition at line 1104 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::cos); }

Field cos (  Field &    a )

Definition at line 764 of file field.cpp. References matvec::Field::map(). {return a.map(std::cos);}

template<class T> Matrix<T> diag (  const Matrix< T > &    a )  [inline]

Definition at line 1099 of file matrix.h. References matvec::Matrix< T >::diag(). { return a.diag(); }

template<class T> Matrix<T> diag (  const Vector< T > &    a )  [inline]

Definition at line 1092 of file matrix.h. References matvec::Vector< T >::size(). Referenced by bound_pd(). { Matrix<T> temp(a.size(),a.size()); for (int i=0; i<a.size(); ++i) temp[i][i] = a[i]; return temp; }

void elemhess (  int    job, double **    mat, const int    n, const int    low, const int    hi, double **    vr, double **    vi )

Definition at line 139 of file nonsymm_eigen.cpp. References check_ptr(). Referenced by nonsymm_eigen(). { //////////////////////////////////////////////////////////////// // Reduce the submatrix in rows and columns low through hi of // real matrix mat to Hessenberg form by elementary similarity // transformations //////////////////////////////////////////////////////////////// int i,j,m; double x,y; int *array = 0; if (hi - low > 1) { if(hi>0){ array = new int [hi]; } else { array = 0; } check_ptr(array); } for (m = low + 1; m < hi; m++) { for (x = 0,i = m,j = m; j <= hi; j++) { if (fabs(mat[j][m-1]) > fabs(x)) { x = mat[j][m-1]; i = j; } } if ((array[m] = i) != m) { for (j = m - 1; j < n; j++) { y = mat[i][j]; mat[i][j] = mat[m][j]; mat[m][j] = y; } for (j = 0; j <= hi; j++) { y = mat[j][i]; mat[j][i] = mat[j][m]; mat[j][m] = y; } } if (x != 0) { for (i = m + 1; i <= hi; i++) { if ((y = mat[i][m-1]) != 0) { y = mat[i][m-1] = y / x; for (j = m; j < n; j++) mat[i][j] -= y * mat[m][j]; for (j = 0; j <= hi; j++) mat[j][m] += y * mat[j][i]; } } } } if (job) { for (i=0; i<n; i++) { for (j=0; j<n; j++) vr[i][j] = 0.0; vi[i][j] = 0.0; vr[i][i] = 1.0; } for (m = hi - 1; m > low; m--) { for (i = m + 1; i <= hi; i++) vr[i][m] = mat[i][m-1]; if ((i = array[m]) != m) { for (j = m; j <= hi; j++) { vr[m][j] = vr[i][j]; vr[i][j] = 0.0; } vr[i][m] = 1.0; } } } if (array) { delete [] array; array=0; } }

template<class T> Matrix<T> erf (  const Matrix< T > &    a )

Definition at line 1115 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(erf); }

Field erf (  Field &    a )

Definition at line 775 of file field.cpp. References matvec::Field::map(). Referenced by Normal_cdf(). {return a.map(::erf);}

template<class T> Matrix<T> erfc (  const Matrix< T > &    a )

Definition at line 1116 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(erfc); }

Field erfc (  Field &    a )

Definition at line 776 of file field.cpp. References matvec::Field::map(). {return a.map(::erfc);}

unsigned matvec::est_nze (  const unsigned    dim, const unsigned    ntrait, const int    ped )

Definition at line 98 of file util.cpp. Referenced by matvec::Model::prepare_data(), matvec::Population::setup_m_ww(), and matvec::Model::vce_emreml_multi_trait(). { // *********************************************** // * get a suggested number of max_nz elements // * (1) calculate max_nz for single trait model // * (2) then multiply it by ntrait*ntrait // ************************************************ if (dim==0) return 0; double onetraitdim = static_cast<double>(dim)/static_cast<double>(ntrait); double nze, halfsize = onetraitdim*(onetraitdim - 1.0)/2.0; if(halfsize < 1000.0) { nze = halfsize; } else if(halfsize < 5000.0) { nze = 0.7 * halfsize; } else if(halfsize < 1.0e+4) { nze = 0.4 * halfsize; } else if(halfsize < 1.0e+5) { nze = 5e-1 * halfsize; } else if(halfsize < 1.0e+6) { nze = 5e-2 * halfsize; } else if(halfsize < 5.0e+6) { nze = 1e-2 * halfsize; } else if(halfsize < 1.0e+7) { nze = 8e-3 * halfsize; } else if(halfsize < 5.0e+7) { nze = 6e-3 * halfsize; } else if(halfsize < 1.0e+8) { nze = 4e-3 * halfsize; } else if(halfsize < 5.0e+8) { nze = 2e-3 * halfsize; } else if(halfsize < 1.0e+9) { nze = 9e-4 * halfsize; } else if(halfsize < 3.0e+9) { nze = 7e-4 * halfsize; } else if(halfsize < 5.0e+9) { nze = 5e-4 * halfsize; } else if(halfsize < 7.0e+9) { nze = 3e-4 * halfsize; } else if(halfsize < 1.0e+10) { nze = 9e-5 * halfsize; } else if(halfsize < 3.0e+10) { nze = 7e-5 * halfsize; } else if(halfsize < 5.0e+10) { nze = 5e-5 * halfsize; } else if(halfsize < 7.0e+10) { nze = 3e-5 * halfsize; } else if(halfsize < 9.0e+10) { nze = 1e-5 * halfsize; } else { nze = 1e-6 * halfsize; } nze += static_cast<double>(pedsize*3); nze = nze*ntrait*(ntrait+1.0)/2.0 + dim; return static_cast<unsigned>(nze); }

template<class T> Matrix<T> exp (  const Matrix< T > &    a )

Definition at line 1112 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::exp); }

Field exp (  Field &    a )

Definition at line 772 of file field.cpp. References matvec::Field::map(). {return a.map(std::exp);}

BG exp (  const BG &    a )

Definition at line 144 of file bg.cpp. References matvec::BG::D, matvec::BG::d, matvec::BG::f, and matvec::Matrix< double >::transpose(). Referenced by Beta_cdf(), betain(), ChiSquare_cdf(), ChiSquare_inv(), factrl(), gamdev(), gammin(), sgamma(), t_cdf(), and t_expected_value(). { BG r; r.f = std::exp(a.f); r.d = r.f*a.d; r.D = r.f*a.D + r.d*a.d.transpose(); return r; }

double matvec::f1dim (  double    x )  [static]

Definition at line 167 of file min_nr_powell.cpp. References funk(), ncom, pcom, and xicom. Referenced by linmin(). { Vector<double> xt(ncom); for (int j=0; j<ncom; j++) xt[j] = pcom[j] + x*xicom[j]; double f = funk(xt); return f; }

double matvec::F_cdf (  const double    f, const double    df1, const double    df2, const double    nc, int &    errcode )

Definition at line 557 of file stat1.cpp. References Beta_cdf(). Referenced by matvec::FDist::cdf(), matvec::GLMM::contrast(), and F_inv(). { if (df1 <= 0.0 || df2 <= 0.0) throw exception(" F_cdf(): bad args: value > 0 expected"); errcode = 0; double retval; if (f <= 0.0) { retval = 0.0; } else { retval = Beta_cdf(df1*f/(df1*f+df2),df1*0.5,df2*0.5,nc,errcode); } return retval; }

double matvec::F_inv (  const double    p, const double    df1, const double    df2, const double    nc )

Definition at line 862 of file stat1.cpp. References F_cdf(). Referenced by matvec::FDist::inv(). { ////////////////////////////////////////////////////////////// // compute percentage point given lower tail area p. // binary iteration algoritm, a better algorithm is needed // // Tianlin Wang, UIUC, Mon Jun 6 14:58:36 CDT 1994 //////////////////////////////////////////////////////////// if (p < 0.0 || p >= 1.0 ) throw exception(" F_inv(): bad arg1, value in (0,1) expected"); if (df1 <= 0.0 || df2 <= 0.0 ) throw exception(" F_inv(): bad arg2 or 3, > 0 expected"); if (p==0.0) return 0.0; double p2,ppf; if (df2 >= 3) { ppf = df2*(df1 + nc)/(df1*(df2 - 2.0)); // mean of F(df1,df2,nc) } else { ppf = 0; } double eps = 1.0e-7; double bot = 0.0; double top = ppf + 26000.0; int errcode = 0; while (top-bot > eps) { ppf = (top + bot)*0.5; p2 = F_cdf(ppf,df1,df2,nc,errcode); if (errcode != 0) break; if (p2 < p ) { bot = ppf; } else if (p2 > p) { top = ppf; } else { break; } } return ppf; }

double matvec::factrl (  const long    n )

Definition at line 315 of file stat1.cpp. References exp(), and gammln(). { static int ntop = 8; static double a[33]={1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0}; int j; if (n > 32) return exp(gammln(static_cast<double>(n+1))); while (ntop < n) { j = ntop++; a[ntop] = a[j]*ntop; } return a[n]; }

int matvec::fieldindex (  const std::string &    fw, const FieldStruct *    field_vec, const unsigned    n )

Definition at line 26 of file fieldstruct.cpp. References matvec::FieldStruct::name(). Referenced by matvec::Model::build_model_term(), matvec::Model::copyfrom(), matvec::Model::covariate(), and matvec::ModelTerm::input(). { int k = -1; for (unsigned i=0; i<n; i++) { if (field_vec[i].name() == fw) { k = i; break; } } return k; }

void flin (  const double    l, const int    j, const Vector< double > &    x, const int    n, const Matrix< double > &    v, Vector< double > &    w, const Matrix< double > &    q01 )  [static]

Definition at line 72 of file min_praxis.cpp. References qd. Referenced by matvec::Minimizer::min_first_dir(). { int i; double qa,qb,qc; if (j != -1) { // linear search for (i=0; i<n; i++) w[i] = x[i] + l *v[i][j]; } else { // search along parabolic space curve qa = l*(l-qd[1])/(qd[0]*(qd[0]+qd[1])); qb = (l+qd[0])*(qd[1]-l)/(qd[0]*qd[1]); qc = l*(l+qd[0])/(qd[1]*(qd[0]+qd[1])); for (i=0; i<n; i++) w[i] = qa*q01[0][i]+qb*x[i]+qc*q01[1][i]; } }

template<class T> Matrix<T> floor (  const Matrix< T > &    a )

Definition at line 1109 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::floor); }

Field floor (  Field &    a )

Definition at line 769 of file field.cpp. References matvec::Field::map(). {return a.map(std::floor);}

double fsign (  double    num, double    sign )

Definition at line 60 of file stat2.cpp. { if ( ( sign>0.0 && num<0.0 ) || ( sign<0.0 && num>0.0 ) ) return ( -num ); else return ( num ); }

double funk (  Vector< double > &    x )  [static]

Definition at line 38 of file min_nr_powell.cpp. References matvec::Model::info(), nfunk, matvec::Model::nterm(), matvec::doubleMatrix::psd(), matvec::Model::residual_var, matvec::Model::restricted_log_likelihood(), matvec::Model::term, and matvec::Model::vec2var(). Referenced by f1dim(), and nr_powell(). { int i,k; unsigned nterms = myM->nterm(); myM->vec2var(x); if (!myM->residual_var.psd()) return 1.0e30; for (i=0; i<nterms; i++) { if (myM->term[i].classi() == 'R' || myM->term[i].classi() == 'P') { k = myM->term[i].prior->var_matrix()->psd(); if (!k) return (1.0e30 + 1.0e10*(i+1)); } } nfunk++; double log_likelihood = myM->restricted_log_likelihood(); if ( (nfunk % 3) == 0) { std::cout << " powell...# of log_likelihood evaluations = " << nfunk << ", " << "and its current value = " << log_likelihood << std::endl; myM->info(std::cout); } return log_likelihood*(-1.0); }

void G_to_gamete (  const unsigned    G_id, unsigned    k1, unsigned    k2, unsigned    size )

Definition at line 41 of file pop_gibbs.cpp. { k1 = G_id/size; k2 = G_id - k1*size; }

double matvec::gamdev (  const double    A, const double    R )

Definition at line 124 of file stat1.cpp. References exp(), log(), ran1(), and sqrt(). { ///////////////////////////////////////////////////////////////////// // Generates random deviates from the standard gamma distribution // whose density is // 1 // --------------- * X^(A-1) * Exp(-X/R) // (R^A)*Gamma(A) // // Chi-Square can be generated as gamdev(df/2,2.0) // // Tianlin Wang at UIUC, May 1,1992 ///////////////////////////////////////////////////////////////////// int K,j; double am,e,s,v1,v2,x,y; K = static_cast<int>(A); if (K < 1) throw exception(" gamdev(): bad arg1, value >= 1 expected "); if (K < 6) { x=1.0; for (j=1;j<=K;j++) x *= ran1(); x = -log(x); } else { do { do { do { v1=2.0*ran1()-1.0; v2=2.0*ran1()-1.0; } while (v1*v1+v2*v2 > 1.0); y=v2/v1; am=K-1; s=sqrt(2.0*am+1.0); x=s*y+am; } while (x <= 0.0); e=(1.0+y*y)*exp(am*log(x/am)-s*y); } while (ran1() > e); } return (x*R); }

double matvec::Gamma_cdf (  const double    g, const double    alfa, const double    theta )

Definition at line 339 of file stat1.cpp. References gammin(). { double retval = 0.0; if (g <= 0.0 ) return retval; if (alfa > 0.0 ) { if (theta > 0.0 ) { retval = gammin(g/theta,alfa); } else { throw exception(" Gamma_cdf(): bad arg3, value > 0 expected"); } } else { throw exception(" Gamma_cdf(): bad arg2: %f, value > 0 expected"); } return retval; }

double matvec::gammin (  const double    x, const double    a )

Definition at line 225 of file stat1.cpp. References exp(), gammln(), log(), Normal_cdf(), and sqrt(). Referenced by ChiSquare_cdf(), ChiSquare_inv(), Gamma_cdf(), and matvec::PoissonDist::sample(). { ////////////////////////////////////////////////////////// // ALGORITHM AS239 Appl. Statist. (1988) 37(3), GAMMAD() // return the incomplete Gamma function for arguments: // x between 0 and +infinite // p positive // REQUIREMENTS: Normal_cdf(), gammln() //////////////////////////////////////////////////////// double zero,one,two,three,nine,tol,oflo,xbig,plimit,elimit; double pn1,pn2,pn3,pn4,pn5,pn6,rn,arg,a,b,c,an; zero = 0.0; one = 1.0; two = 2.0; three = 3.0; nine = 9.0; tol = 1.0e-14; oflo = 1.0e+37; xbig = 1.0e+8; plimit = 1000.0; elimit = -88.0; double retval = zero; if (x < zero || p <= zero) throw exception(" gammin(): bad args"); if (x == zero) return retval; if (p > plimit) { // use a normal approximation pn1 = three*sqrt(p)*(pow(x/p,one/three) + one/(nine*p) - one); return Normal_cdf(pn1); } if (x > xbig) return one; // x is extremely large compared to p if (x <= one || x < p) { // use Pearson's series expansion arg = p*log(x) - x - gammln(p + one); c = one; retval = one; a = p; do { a++; c *= x/a; retval += c; } while (c > tol); arg += log(retval); retval = zero; if (arg >= elimit) retval = exp(arg); } else { arg = p*log(x) - x - gammln(p); a = one - p; b = a + x + one; c = zero; pn1 = one; pn2 = x; pn3 = x + one; pn4 = x*b; retval = pn3/pn4; for(;;) { a += one; b += two; c += one; an = a*c; pn5 = b*pn3 - an*pn1; pn6 = b*pn4 - an*pn2; if (fabs(pn6) > zero) { rn = pn5/pn6; if (fabs(retval - rn) <= tol*(rn < one ? rn: one)) break; retval = rn; } pn1 = pn3; pn2 = pn4; pn3 = pn5; pn4 = pn6; if (fabs(pn5) >= oflo) { // re-scale terms, if terms are large pn1 /= oflo; pn2 /= oflo; pn3 /= oflo; pn4 /= oflo; } } arg += log(retval); retval = one; if (arg >= elimit) retval = one - exp(arg); } return retval; }

double matvec::gammln (  const double    xx )

Definition at line 208 of file stat1.cpp. References log(). Referenced by Beta_cdf(), BinCoef(), matvec::BinomialDist::cdf(), matvec::GammaDist::cdf(), ChiSquare_cdf(), ChiSquare_inv(), factrl(), gammin(), matvec::NegativeBinomialDist::mgf(), matvec::BinomialDist::sample(), t_cdf(), and t_expected_value(). { if (xx <= 0.0) throw exception(" gammln(): bad arg, value > 0 expected"); double x,tmp,ser; static double cof[6]={76.18009173,-86.50532033,24.01409822, -1.231739516,0.120858003e-2,-0.536382e-5}; x = xx-1.0; tmp = x + 5.5; tmp -= (x+0.5)*log(tmp); ser=1.0; for (int j=0; j<=5; j++) { x += 1.0; ser += cof[j]/x; } return -tmp+log(ser*2.50662827465); }

double matvec::gasdev (  void    )

Definition at line 103 of file stat1.cpp. References log(), ran1(), and sqrt(). { static int iset=0; static double gset; double fac,r,v1,v2; if (iset == 0) { do { v1=2.0*ran1()-1.0; v2=2.0*ran1()-1.0; r=v1*v1+v2*v2; } while (r >= 1.0 || r == 0.0); fac=sqrt(-2.0*log(r)/r); gset=v1*fac; iset=1; return v2*fac; } else { iset=0; return gset; } }

double matvec::genbet (  const double    aa, const double    bb )

Definition at line 1135 of file stat2.cpp. References genbet(), min, and ranf(). Referenced by genbet(), and matvec::BetaDist::sample(). :317-322 (1978) // (Algorithms BB and BC) // *********************************************************************** { #define expmax 89.0 #define infnty 1.0E38 static double olda = -1.0; static double oldb = -1.0; static double genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z; static long qsame; qsame = olda == aa && oldb == bb; if (qsame) goto S20; if (aa <= 0.0 || bb <= 0.0) throw exception("genbet(): bad args"); olda = aa; oldb = bb; S20: if(!(std::min(aa,bb) > 1.0)) goto S100; // Alborithm BB Initialize if(qsame) goto S30; a = std::min(aa,bb); b = std::max(aa,bb); alpha = a+b; beta = sqrt((alpha-2.0)/(2.0*a*b-alpha)); gamma = a+1.0/beta; S30: S40: u1 = ranf(); // Step 1 u2 = ranf(); v = beta*log(u1/(1.0-u1)); if(!(v > expmax)) goto S50; w = infnty; goto S60; S50: w = a*exp(v); S60: z = pow(u1,2.0)*u2; r = gamma*v-1.3862944; s = a+r-w; // Step 2 if(s+2.609438 >= 5.0*z) goto S70; // Step 3 t = log(z); if(s > t) goto S70; // Step 4 if(r+alpha*log(alpha/(b+w)) < t) goto S40; S70: // Step 5 if(!(aa == a)) goto S80; genbet = w/(b+w); goto S90; S80: genbet = b/(b+w); S90: goto S230; S100: // Algorithm BC Initialize if(qsame) goto S110; a = std::max(aa,bb); b = std::min(aa,bb); alpha = a+b; beta = 1.0/b; delta = 1.0+a-b; k1 = delta*(1.38889e-2+4.16667e-2*b)/(a*beta-0.777778); k2 = 0.25+(0.5+0.25/delta)*b; S110: S120: u1 = ranf(); // Step 1 u2 = ranf(); if(u1 >= 0.5) goto S130; //Step 2 y = u1*u2; z = u1*y; if(0.25*u2+z-y >= k1) goto S120; goto S170; S130: // Step 3 z = pow(u1,2.0)*u2; if(!(z <= 0.25)) goto S160; v = beta*log(u1/(1.0-u1)); if(!(v > expmax)) goto S140; w = infnty; goto S150; S140: w = a*exp(v); S150: goto S200; S160: if(z >= k2) goto S120; S170: // Step 4 Step 5 v = beta*log(u1/(1.0-u1)); if(!(v > expmax)) goto S180; w = infnty; goto S190; S180: w = a*exp(v); S190: if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120; S200: // Step 6 if(!(a == aa)) goto S210; genbet = w/(b+w); goto S220; S210: genbet = b/(b+w); S230: S220: return genbet; #undef expmax #undef infnty }

double matvec::genchi (  const int    df )

Definition at line 487 of file stat2.cpp. References genchi(), and sgamma(). Referenced by genchi(), genf(), gennch(), gennf(), matvec::tDist::sample(), matvec::ChiSquareDist::sample(), and matvec::Model::vce_gibbs_sampler(). { if (df <= 0) throw exception(" genchi(): invalid arg"); return (2.0*sgamma(static_cast<double>(df)/2.0)); }

double matvec::genexp (  const double    av )

Definition at line 508 of file stat2.cpp. References genexp(), and sexpo(). Referenced by genexp(). : // Ahrens, J.H. and Dieter, U. // Computer Methods for Sampling From the // Exponential and Normal Distributions. // Comm. ACM, 15,10 (Oct. 1972), 873 - 882. //************************************************************ { return ( sexpo()*av ); }

double matvec::genf (  const int    dfn, const int    dfd )

Definition at line 534 of file stat2.cpp. References genchi(), genf(), and warning(). Referenced by genf(), and matvec::FDist::sample(). { if ( dfn <= 0.0 || dfd <= 0.0) throw exception(" genf(): invalid args"); static double genfval = 0.0; static double xden,xnum; xnum = genchi(dfn)/dfn; ////// GENF = ( GENCHI(DFN)/DFN )/( GENCHI(DFD)/DFD ) /////// xden = genchi(dfd)/dfd; if ( xden <= 9.999999999998e-39*xnum) { warning("GENF - generated numbers would cause overflow," " df1 = %16.6E, df2 = %16.6E",xnum,xden); genfval = 1.0E38; } else { genfval = xnum/xden; } return genfval; }

double matvec::gengam (  const double    a, const double    r )

Definition at line 451 of file stat2.cpp. References gengam(), and sgamma(). Referenced by gengam(). : // (Case A >= 1.0) // Ahrens, J.H. and Dieter, U. // Generating Gamma Variates by a // Modified Rejection Technique. // Comm. ACM, 25,1 (Jan. 1982), 47 - 54. // Algorithm GD // (Case 0.0 <= A <= 1.0) // Ahrens, J.H. and Dieter, U. // Computer Methods for Sampling from Gamma, // Beta, Poisson and Binomial Distributions. // Computing, 12 (1974), 223-246/ // Adapted algorithm GS. //***************************************************************** { return ( sgamma(a)*r ); }

double matvec::gennch (  const int    df, const double    xnonc )

Definition at line 570 of file stat2.cpp. References genchi(), gennch(), gennor(), and snorm(). Referenced by gennch(), and gennf(). { if ( df > 1 && xnonc >= 0.0) { return ( genchi(df-1)+pow(gennor(sqrt(xnonc),1.0),2.0)); } else if (df == 1 && xnonc >= 0.0) { double x = snorm() + xnonc; return x*x; } else { throw exception(" gennch(): invalid arg"); } }

double matvec::gennf (  const int    dfn, const int    dfd, const double    xnonc )

Definition at line 602 of file stat2.cpp. References genchi(), gennch(), gennf(), and warning(). Referenced by gennf(), and matvec::FDist::sample(). { if (dfn <= 1 || dfd <= 0 || xnonc < 0.0) throw exception(" gennf(): invalid args"); static double gennf,xden,xnum; xnum = gennch(dfn,xnonc)/dfn; ///// GENNF = ( GENNCH(DFN,XNONC)/DFN )/( GENCHI(DFD)/DFD ) ///// xden = genchi(dfd)/dfd; if ( xden <= 9.999999999998e-39*xnum ) { warning("gennf(): overflow, df1= %16.6E, df2= %16.6E",xnum,xden); gennf = 1.0E38; } else { gennf = xnum/xden; } return gennf; }

double matvec::gennor (  const double    av, const double    sd )

Definition at line 428 of file stat2.cpp. References gennor(), and snorm(). Referenced by gennch(), and gennor(). : // Ahrens, J.H. and Dieter, U. // Extensions of Forsythe's Method for Random // Sampling from the Normal Distribution. // Math. Comput., 27,124 (Oct. 1973), 927 - 937. //******************************************************************** { return ( sd*snorm()+av ); }

double matvec::genunf (  const double    low, const double    high )

Definition at line 638 of file stat2.cpp. References genunf(), and ranf(). Referenced by genunf(), and matvec::UniformDist::sample(). { if (low > high ) throw exception("genunf(): invalid arg"); return ( low+(high-low)*ranf() ); }

void matvec::getlambda (  double **    lambda, const int    n )

Definition at line 37 of file pedigree.cpp. Referenced by matvec::Model::add_Ag(), matvec::GLMM::add_Ag(), matvec::Model::add_Ag_diag(), matvec::GLMM::add_Ag_old(), matvec::GLMM::add_AgSand(), matvec::Model::add_G_1_single_trait(), matvec::Population::add_Ga_inv(), matvec::Pedigree::ainv(), and matvec::GLMM::ainv(). { if (lambda) { if (n == 3 ) { lambda[0][0] = 1.0; lambda[0][1] = -0.5; lambda[0][2] = -0.5; lambda[1][0] = -0.5; lambda[1][1] = 0.25; lambda[1][2] = 0.25; lambda[2][0] = -0.5; lambda[2][1] = 0.25; lambda[2][2] = 0.25; } else { throw exception("getlambda(a,n),n must be 3"); } } else { throw exception("getlambda(l,n): l is NULL"); } }

unsigned matvec::ginverse1 (  BG **    a, const unsigned    n, BG &    lgdet, int    mode, const double    tol )

Definition at line 112 of file ginverse1.cpp. References ginverse1(). { unsigned i,j,k,irank = n; BG sum; /////////////////////////////////////////////////// // column-by-column, we are doing // L (lower triangular) in place such that L*L'=A /////////////////////////////////////////////////// lgdet = 0.0; for (j=0; j<n; ++j) { // for column j, L[j][j] is calculated first sum = a[j][j]; for (k=0; k<j; ++k) sum -= a[j][k]*a[j][k]; // for other elements L[i,j] in jth column, i=j+1 to n if (sum <= -tol) { throw exception("ginverse1(): matrix is non-psd"); } if (fabs(sum) < tol) { irank--; for (i=j; i<n; i++) a[i][j]=0.0; } else { lgdet += log(sum); a[j][j]=sqrt(sum); for (i=j+1; i<n; i++) { sum = a[i][j]; for (k=0; k<j; ++k) sum -= a[i][k]*a[j][k]; a[i][j]=sum/a[j][j]; } } } if (mode != 1) { if (mode ==2) { for(i=0;i<n-1;i++) { //memset(&(a[i][i+1]),'\0',sizeof(BG)*(n-i-1)); for(j=i+1;j<n;j++){ a[i][j] = 0.0; } } } return irank; } ////////////////////////////////////////// // column-by-column, we are doing // inv(L) in place (lower triangl) ////////////////////////////////////////// for (j=0; j<n; j++) { if (fabs(a[j][j]) < tol) { for (i=j; i<n; i++) a[i][j]=0.0; } else { a[j][j] = 1./a[j][j]; for (i=j+1; i<n; i++) { if (a[i][i] <= -tol) throw exception("ginverse1(): matrix is non-psd"); if (fabs(a[i][i]) < tol) { a[i][j]=0.0; } else { sum=0.0; for (k=j; k<i; k++) sum -= a[i][k]*a[k][j]; a[i][j] = sum/a[i][i]; } // end of if-else block } // end of for-loop } // end of if-else block } // end of for-loop ////////////////////////////////////// // column-by-column, we are doing // inv(L')*inv(L)= inv(a) in place ////////////////////////////////////// for (j=0; j<n; j++) { for (i=j; i<n; i++) { for (sum=0.0,k=i; k<n; k++) sum += a[k][i]*a[k][j]; a[i][j]=sum; a[j][i]=sum; } } return irank; }

unsigned matvec::ginverse1 (  double **    a, const unsigned    n, double &    lgdet, int    mode, const double    tol )

Definition at line 38 of file ginverse1.cpp. Referenced by matvec::Model::add_Ag(), matvec::GLMM::add_Ag(), matvec::Model::add_Ag_diag(), matvec::GLMM::add_Ag_old(), matvec::GLMM::add_AgSand(), matvec::Model::add_Ig(), matvec::Model::add_Ig_diag(), matvec::GLMM::add_IgSand(), matvec::doubleMatrix::ginv1(), matvec::BGMatrix::ginv1(), ginverse1(), matvec::Model::inverse_residual_var(), matvec::doubleMatrix::logdet(), matvec::BGMatrix::logdet(), matvec::doubleMatrix::sqrtm(), and matvec::BGMatrix::sqrtm(). { unsigned i,j,k,irank = n; double sum; /////////////////////////////////////////////////// // column-by-column, we are doing // L (lower triangular) in place such that L*L'=A /////////////////////////////////////////////////// lgdet = 0.0; for (j=0; j<n; ++j) { // for column j, L[j][j] is calculated first sum = a[j][j]; for (k=0; k<j; ++k) sum -= a[j][k]*a[j][k]; // for other elements L[i,j] in jth column, i=j+1 to n if (sum <= -tol) { throw exception("ginverse1(): matrix is non-psd"); } if (std::fabs(sum) < tol) { irank--; for (i=j; i<n; i++) a[i][j]=0.0; } else { lgdet += std::log(sum); a[j][j]=std::sqrt(sum); for (i=j+1; i<n; i++) { sum = a[i][j]; for (k=0; k<j; ++k) sum -= a[i][k]*a[j][k]; a[i][j]=sum/a[j][j]; } } } if (mode != 1) { if (mode ==2) { for(i=0;i<n-1;i++) memset(&(a[i][i+1]),'\0',sizeof(double)*(n-i-1)); } return irank; } ////////////////////////////////////////// // column-by-column, we are doing // inv(L) in place (lower triangl) ////////////////////////////////////////// for (j=0; j<n; j++) { if (std::fabs(a[j][j]) < tol) { for (i=j; i<n; i++) a[i][j]=0.0; } else { a[j][j] = 1./a[j][j]; for (i=j+1; i<n; i++) { if (a[i][i] <= -tol) throw exception("ginverse1(): matrix is non-psd"); if (std::fabs(a[i][i]) < tol) { a[i][j]=0.0; } else { sum=0.0; for (k=j; k<i; k++) sum -= a[i][k]*a[k][j]; a[i][j] = sum/a[i][i]; } // end of if-else block } // end of for-loop } // end of if-else block } // end of for-loop ////////////////////////////////////// // column-by-column, we are doing // inv(L')*inv(L)= inv(a) in place ////////////////////////////////////// for (j=0; j<n; j++) { for (i=j; i<n; i++) { for (sum=0.0,k=i; k<n; k++) sum += a[k][i]*a[k][j]; a[i][j]=sum; a[j][i]=sum; } } return irank; }

void ginverse2 (  double *    a, const unsigned    n, const double    tol )  [static]

Definition at line 40 of file model_blup.cpp. References lidx(). Referenced by matvec::Model::blup_pccg(). { unsigned i,j,k,t,ii,jj; double sum,aii,ajj; /////////////////////////////////////////////////// // column-by-column, we are doing // L (lower triangular) in place such that L*L'=A /////////////////////////////////////////////////// for (j=0; j<n; j++) { // for column j, L[j][j] is calculated first jj = lidx(j,j,n); sum = ajj = a[jj]; for (k=0; k<j; ++k) { aii = a[lidx(j,k,n)]; sum -= aii*aii; } // for other elements L[i,j] in jth column, i=j+1 to n if (sum <= -tol) { throw exception("ginverse2(): matrix is non-psd"); } if (fabs(sum) < tol) { for (i=j; i<n; i++) { a[jj++]=0.0; } } else { a[jj]= ajj = std::sqrt(sum); for (i=j+1; i<n; i++) { sum = a[lidx(i,j,n)]; for (k=0; k<j; ++k) { sum -= a[lidx(i,k,n)]*a[lidx(j,k,n)]; } a[lidx(i,j,n)]=sum/ajj; } } } ////////////////////////////////////////// // column-by-column, we are doing // inv(L) in place (lower triangl) ////////////////////////////////////////// for (j=0; j<n; j++) { jj = lidx(j,j,n); if (a[jj] < tol) { for (i=j; i<n; i++) { a[jj++]=0.0; } } else { a[jj] = 1./a[jj]; for (i=j+1; i<n; i++) { aii = a[lidx(i,i,n)]; jj++; if (aii <= -tol) { throw exception("ginverse2(): matrix is non-psd"); } else if (fabs(aii) < tol) { a[jj]=0.0; } else { for (sum=0.0,k=j; k<i; ++k) sum -= a[lidx(i,k,n)]*a[lidx(k,j,n)]; a[jj] = sum/aii; } // end of if-else block } // end of for-loop } // end of if-else block } // end of for-loop ////////////////////////////////////// // column-by-column, we are doing // inv(L')*inv(L)= inv(a) in place ////////////////////////////////////// for (j=0; j<n; j++) { jj = lidx(j,j,n); for (i=j; i<n; i++) { ii = lidx(i,i,n); k = lidx(i,j,n); for (sum=0.0,t=i; t<n; ++t) sum += a[ii++]*a[k++]; a[jj++]=sum; } } return; }

void matvec::gs_ioc (  int *    ia, int *    ja, double *    a, const double *    rhs, double *    sol, const int    n, const double    relax, const double    stopval, const int    mxiter, const double    tol, unsigned &    cov )

Definition at line 55 of file sparsematrix.cpp. { int j,row,col; double oldsol,newsol,local,diag, cmax,cval; Vector<double> adj_rhs(n); int niter = 0; cov = 0; std::cout << std::endl; do { // now iteration begins cmax = 0.0; for (row=1; row<=n; row++) adj_rhs[row-1] = rhs[row]; for (row=1; row<=n; row++) { oldsol = sol[row]; local = 0.0; diag = 0.0; // subtract all the other solutions*count: for (j=ia[row]; j<=ia[row+1]-1; j++) { col = ja[j]; if (col != row) { local += a[j]*sol[col]; } else { diag = a[j]; } } if (diag > tol) { newsol = (adj_rhs[row - 1] - local)/diag; // new solution for GS cval = (newsol - oldsol)*relax; // change vs old solution if (fabs(cval) > cmax) cmax=fabs(cval); newsol = sol[row] = oldsol + cval; // new solution for SOR newsol = sol[row]; for (j = ia[row]; j<=ia[row+1]-1; j++) { col = ja[j]; adj_rhs[col - 1] -= a[j]*newsol; // adjust the adj_rhs } } else if (diag > -tol) { throw exception(" zero-diagonal found in gs_ioc.c"); } else { throw exception(" in gs_ioc(), matrix is not psd"); } } niter += 1; if ( (niter % 10) ==0) { std::cout << " IOC: # of iter = " << niter << ", max_change = " << cmax << std::endl; } } while (niter <= mxiter && cmax > stopval); if (cmax <= stopval) cov=1; }

void matvec::gs_ioc (  int *    ia, int *    ja, BG *    a, const BG *    rhs, BG *    sol, const int    n, const double    relax, const double    stopval, const int    mxiter, const double    tol, unsigned &    cov )

Definition at line 55 of file sparsebgmatrix.cpp. References matvec::BG::f. Referenced by matvec::SparseMatrix::solve(), and matvec::SparseBGMatrix::solve(). { int j,row,col; BG oldsol,newsol,local,diag, cmax,cval; Vector<BG> adj_rhs(n); int niter = 0; cov = 0; std::cout << std::endl; do { // now iteration begins cmax = 0.0; for (row=1; row<=n; row++) adj_rhs[row-1] = rhs[row]; for (row=1; row<=n; row++) { oldsol = sol[row]; local = 0.0; diag = 0.0; // subtract all the other solutions*count: for (j=ia[row]; j<=ia[row+1]-1; j++) { col = ja[j]; if (col != row) { local += a[j]*sol[col]; } else { diag = a[j]; } } if (diag > tol) { newsol = (adj_rhs[row - 1] - local)/diag; // new solution for GS cval = (newsol - oldsol)*relax; // change vs old solution if (fabs(cval) > cmax) cmax=fabs(cval); newsol = sol[row] = oldsol + cval; // new solution for SOR newsol = sol[row]; for (j = ia[row]; j<=ia[row+1]-1; j++) { col = ja[j]; adj_rhs[col - 1] -= a[j]*newsol; // adjust the adj_rhs } } else if (diag > -tol) { throw exception(" zero-diagonal found in gs_ioc.c"); } else { throw exception(" in gs_ioc(), matrix is not psd"); } } niter += 1; if ( (niter % 10) ==0) { std::cout << " IOC: # of iter = " << niter << ", max_change = " << cmax.f << std::endl; } } while (niter <= mxiter && cmax > stopval); if (cmax <= stopval) cov=1; }

template<class T> Matrix<T> hadjoin (  const Vector< T > &    a, const Matrix< T > &    b )  [inline]

Definition at line 1048 of file matrix.h. References matvec::Matrix< T >::num_cols(), matvec::Matrix< T >::num_rows(), and matvec::Vector< T >::size(). { if (a.size() != b.num_rows()) throw exception("Matrix<T>::hadjoin(): size not conformable"); Matrix<T> temp(a.size(),1 + b.num_cols()); for (int i=0; i<a.size(); ++i) { temp[i][0] = a[i]; memcpy(&(temp[i][1]),b[i],sizeof(T)*b.num_cols()); } return temp; }

template<class T> Matrix<T> hadjoin (  const Matrix< T > &    a, const Vector< T > &    b )  [inline]

Definition at line 1037 of file matrix.h. References matvec::Matrix< T >::num_cols(), matvec::Matrix< T >::num_rows(), and matvec::Vector< T >::size(). { if (a.num_rows() != b.size()) throw exception("Matrix<T>::hadjoin(): size not conformable"); Matrix<T> temp(a.num_rows(),a.num_cols() + 1); for (int i=0; i<a.num_rows(); ++i) { memcpy(temp[i],a[i],sizeof(T)*a.num_cols()); temp[i][a.num_cols()] = b[i]; } return temp; }

template<class T> Matrix<T> hadjoin (  const Matrix< T > &    a, const Matrix< T > &    b )  [inline]

Definition at line 1026 of file matrix.h. References matvec::Matrix< T >::num_cols(), and matvec::Matrix< T >::num_rows(). Referenced by matvec::GLMM::AI_REML(). { if (a.num_rows() != b.num_rows()) throw exception("Matrix<T>::hadjoin(): size not conformable"); Matrix<T> temp(a.num_rows(),a.num_cols() + b.num_cols()); for (int i=0; i<a.num_rows(); ++i) { memcpy(temp[i],a[i],sizeof(T)*a.num_cols()); memcpy(&(temp[i][a.num_cols()]),b[i],sizeof(T)*b.num_cols()); } return temp; }

void matvec::hsort_ija (  unsigned    n, int *    ia, int *    ja, BG *    a )  [static]

this is a Heap sort, sorting (ia,ja) ascendingly. ia(1...n), ja(1...n) ,a(1..n) where n is the actual number of non-zeros in A Definition at line 1231 of file sparsebgmatrix.cpp. { unsigned iia, jja,i,j,l,ir; BG aa; if (n == 1) return; l = (n >> 1)+1; /* l = n/2+1; */ ir=n; for (;;) { if (l > 1) { l = l-1; iia = ia[l]; jja = ja[l]; aa = a[l]; } else { iia = ia[ir]; jja = ja[ir]; aa = a[ir]; ia[ir] = ia[1]; ja[ir] = ja[1]; a[ir] = a[1]; ir--; if(ir == 1) { ia[1] = iia; ja[1] = jja; a[1] = aa; return; } } i=l; j=l+l; while ( j <= ir) { if (j < ir) { if (ia[j]<ia[j+1] || (ia[j]==ia[j+1] && ja[j] < ja[j+1])) j++; } if (iia < ia[j] || (iia == ia[j] && jja < ja[j])) { ia[i] = ia[j]; ja[i] = ja[j]; a[i] = a[j]; i=j; j=j+j; } else { j=ir+1; } } ia[i] = iia; ja[i] = jja; a[i] = aa; } }

void hsort_ija (  unsigned    n, int *    ia, int *    ja, double *    a )  [static]

this is a Heap sort, sorting (ia,ja) ascendingly. ia(1...n), ja(1...n) ,a(1..n) where n is the actual number of non-zeros in A Definition at line 1642 of file glmm2.cpp. Referenced by matvec::GLMM::build_hInv(), matvec::SparseMatrix::close(), and matvec::SparseBGMatrix::close(). { unsigned iia, jja,i,j,l,ir; double aa; if (n == 1) return; l = (n >> 1)+1; /* l = n/2+1; */ ir=n; for (;;) { if (l > 1) { l = l-1; iia = ia[l]; jja = ja[l]; aa = a[l]; } else { iia = ia[ir]; jja = ja[ir]; aa = a[ir]; ia[ir] = ia[1]; ja[ir] = ja[1]; a[ir] = a[1]; ir--; if(ir == 1) { ia[1] = iia; ja[1] = jja; a[1] = aa; return; } } i=l; j=l+l; while ( j <= ir) { if (j < ir) { if (ia[j]<ia[j+1] || (ia[j]==ia[j+1] && ja[j] < ja[j+1])) j++; } if (iia < ia[j] || (iia == ia[j] && jja < ja[j])) { ia[i] = ia[j]; ja[i] = ja[j]; a[i] = a[j]; i=j; j=j+j; } else { j=ir+1; } } ia[i] = iia; ja[i] = jja; a[i] = aa; } }

long matvec::ignbin (  const long    n, const double    pp )

Definition at line 653 of file stat2.cpp. References ignbin(), min, and ranf(). Referenced by ignbin(), matvec::BinomialDist::sample(), and matvec::Pedigree::sample(). : // Kachitvichyanukul, V. and Schmeiser, B. W. // Binomial Random Variate Generation. // Communications of the ACM, 31, 2 // (February, 1988) 216. // // SUBROUTINE BTPEC(N,PP,ISEED,JX) // BINOMIAL RANDOM VARIATE GENERATOR // MEAN .LT. 30 -- INVERSE CDF // MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA // FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE // (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE // THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS. // BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL. // BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE // RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE // USABLE ALGORITHM. // REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER, // "BINOMIAL RANDOM VARIATE GENERATION," // COMMUNICATIONS OF THE ACM, FORTHCOMING // WRITTEN: SEPTEMBER 1980. // LAST REVISED: MAY 1985, JULY 1987 // REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER // GENERATOR // ARGUMENTS // N : NUMBER OF BERNOULLI TRIALS (INPUT) // PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT) // ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT) // JX: RANDOMLY GENERATED OBSERVATION (OUTPUT) // VARIABLES // PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC // NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC // XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC // P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC // FFM: TEMPORARY VARIABLE EQUAL TO XNP + P // M: INTEGER VALUE OF THE CURRENT MODE // FM: FLOATING POINT VALUE OF THE CURRENT MODE // XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS // P1: AREA OF THE TRIANGLE // C: HEIGHT OF THE PARALLELOGRAMS // XM: CENTER OF THE TRIANGLE // XL: LEFT END OF THE TRIANGLE // XR: RIGHT END OF THE TRIANGLE // AL: TEMPORARY VARIABLE // XLL: RATE FOR THE LEFT EXPONENTIAL TAIL // XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL // P2: AREA OF THE PARALLELOGRAMS // P3: AREA OF THE LEFT EXPONENTIAL TAIL // P4: AREA OF THE RIGHT EXPONENTIAL TAIL // U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE // FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE // FROM THE REGION // V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE // (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR // REJECT THE CANDIDATE VALUE // IX: INTEGER CANDIDATE VALUE // X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC // AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC // K: ABSOLUTE VALUE OF (IX-M) // F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE // ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL // ALSO USED IN THE INVERSE TRANSFORMATION // R: THE RATIO P/Q // G: CONSTANT USED IN CALCULATION OF PROBABILITY // MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION // OF F WHEN IX IS GREATER THAN M // IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT // CALCULATION OF F WHEN IX IS LESS THAN M // I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE // AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND // YNORM: LOGARITHM OF NORMAL BOUND // ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V // X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE // USED IN THE FINAL ACCEPT/REJECT TEST // QN: PROBABILITY OF NO SUCCESS IN N TRIALS // REMARK // IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD // SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME // COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS // ARE NOT INVOLVED. // ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE // GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE // TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR // *DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY* //************************************************************************ { if (pp < 0.0 || pp>1.0) throw exception("ignbin(): bad arg2, value between (0,1) expected"); static double psave = -1.0; static long nsave = -1; static long ignbin,i,ix,ix1,k,m,mp,T1; static double al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,q,qn,r; static double u,v,w,w2,x,x1,x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2; if(pp != psave) goto S10; if(n != nsave) goto S20; if(xnp < 30.0) goto S150; goto S30; S10: /////// SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE ////// psave = pp; p = std::min(psave,1.0-psave); q = 1.0-p; S20: xnp = n*p; nsave = n; if(xnp < 30.0) goto S140; ffm = xnp+p; m = long(ffm); fm = m; xnpq = xnp*q; p1 = (long) (2.195*sqrt(xnpq)-4.6*q)+0.5; xm = fm+0.5; xl = xm-p1; xr = xm+p1; c = 0.134+20.5/(15.3+fm); al = (ffm-xl)/(ffm-xl*p); xll = al*(1.0+0.5*al); al = (xr-ffm)/(xr*q); xlr = al*(1.0+0.5*al); p2 = p1*(1.0+c+c); p3 = p2+c/xll; p4 = p3+c/xlr; S30: //////// GENERATE VARIATE //////// u = ranf()*p4; v = ranf(); /////// TRIANGULAR REGION /////// if(u > p1) goto S40; ix = long(xm-p1*v+u); goto S170; S40: ///// PARALLELOGRAM REGION ///// if(u > p2) goto S50; x = xl+(u-p1)/c; v = v*c+1.0-fabs(xm-x)/p1; if(v > 1.0 || v <= 0.0) goto S30; ix = long(x); goto S70; S50: ////// LEFT TAIL /////// if(u > p3) goto S60; ix = long(xl+log(v)/xll); if(ix < 0) goto S30; v *= ((u-p2)*xll); goto S70; S60: //// RIGHT TAIL //// ix = long(xr-log(v)/xlr); if(ix > n) goto S30; v *= ((u-p3)*xlr); S70: //// DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST //// k = abs((int)(ix-m)); if(k > 20 && k < xnpq/2-1) goto S130; /////// EXPLICIT EVALUATION ////// f = 1.0; r = p/q; g = (n+1)*r; T1 = m-ix; if(T1 < 0) { goto S80; } else if(T1 == 0) { goto S120; } else { goto S100; } S80: mp = m+1; for(i=mp; i<=ix; i++) f *= (g/i-r); goto S120; S100: ix1 = ix+1; for(i=ix1; i<=m; i++) f /= (g/i-r); S120: if(v <= f) goto S170; goto S30; S130: ///// SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X)) //// amaxp = k/xnpq*((k*(k/3.0+0.625)+0.1666666666666)/xnpq+0.5); ynorm = -(k*k/(2.0*xnpq)); alv = log(v); if(alv < ynorm-amaxp) goto S170; if(alv > ynorm+amaxp) goto S30; /////////////////////////////////////////////////// // STIRLING'S FORMULA TO MACHINE ACCURACY FOR // THE FINAL ACCEPTANCE/REJECTION TEST //////////////////////////////////////////////////// x1 = ix+1.0; f1 = fm+1.0; z = n+1.0-fm; w = n-ix+1.0; z2 = z*z; x2 = x1*x1; f2 = f1*f1; w2 = w*w; if (alv <= xm*log(f1/x1)+(n-m+0.5)*log(z/w)+(ix-m)*log(w*p/(x1*q))+(13860.0- (462.0-(132.0-(99.0-140.0/f2)/f2)/f2)/f2)/f1/166320.0+(13860.0-(462.0- (132.0-(99.0-140.0/z2)/z2)/z2)/z2)/z/166320.0+(13860.0-(462.0-(132.0- (99.0-140.0/x2)/x2)/x2)/x2)/x1/166320.0+(13860.0-(462.0-(132.0-(99.0 -140.0/w2)/w2)/w2)/w2)/w/166320.0) goto S170; goto S30; S140: //// INVERSE CDF LOGIC FOR MEAN LESS THAN 30 //// qn = pow(q,(double)n); r = p/q; g = r*(n+1); S150: ix = 0; f = qn; u = ranf(); S160: if(u < f) goto S170; if(ix > 110) goto S150; u -= f; ix += 1; f *= (g/ix-r); goto S160; S170: if(psave > 0.5) ix = n-ix; ignbin = ix; return ignbin; }

long matvec::ignpoi (  const double    mu )

Definition at line 894 of file stat2.cpp. References ignpoi(), min, and ranf(). Referenced by ignpoi(), and matvec::PoissonDist::sample(). { ///////////////////////////////////////////////////////////////////////// // all labels stattements have been replaced with non-labeled statements // T Wang, 4/27/94 ///////////////////////////////////////////////////////////////////////// static double a0 = -0.5; static double a1 = 0.3333333; static double a2 = -0.2500068; static double a3 = 0.2000118; static double a4 = -0.1661269; static double a5 = 0.1421878; static double a6 = -0.1384794; static double a7 = 0.125006; static double muold = 0.0; static double muprev = 0.0; static double fact[10] = { 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0 }; static long ignpoi,j,k,kflag,l,m; static double b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,fk,fx,fy,g; static double omega,p,p0,px,py,q,s,t,u,v,x,xx,pp[35]; int wflag; // added by T. Wang if (mu != muprev) { if (mu < 10.0) { // C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY) muprev = 0.0; if (mu != muold) { muold = mu; m = std::max(1L,(long) (mu)); l = 0; p = exp(-mu); q = p0 = p; } // STEP U. UNIFORM SAMPLE FOR INVERSION METHOD for (;;) { u = ranf(); ignpoi = 0; if (u <= p0) return ignpoi; ////////////////////////////////////////////////////////// // STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE // PP-TABLE OF CUMULATIVE POISSON PROBABILITIES // (0.458=PP(9) FOR MU=10) ////////////////////////////////////////////////////////// if (l == 0) { //////////////////////////////////////////////////////// // STEP C. CREATION OF NEW POISSON PROBABILITIES P // AND THEIR CUMULATIVES Q=PP(K) ///////////////////////////////////////////////////////// l += 1; for (k=l; k<=35; k++) { p = p*mu/(double)k; q += p; *(pp+k-1) = q; if(u <= q) {l = ignpoi = k; return ignpoi;} } l = 35; } else { j = 1; if (u > 0.458) j = std::min(l,m); for (k=j; k<=l; k++) { if (u <= *(pp+k-1)) {ignpoi = k; return ignpoi;} } if (l != 35) { l += 1; for (k=l; k<=35; k++) { p = p*mu/(double)k; q += p; *(pp+k-1) = q; if(u <= q) {l = ignpoi = k; return ignpoi;} } l = 35; } } } } else { //// C A S E A. (RECALCULATION OF S,D,L IF MU HAS CHANGED) //// muprev = mu; s = sqrt(mu); d = 6.0*mu*mu; ///////////////////////////////////////////////////////////// // THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL // PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484) // IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 . /////////////////////////////////////////////////////////////// l = (long) (mu-1.1484); } } //// STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE /// g = mu+s*snorm(); if (g >= 0.0) { ignpoi = (long) (g); // STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH // if (ignpoi >= l) return ignpoi; // STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U // fk = (double)ignpoi; difmuk = mu-fk; u = ranf(); if (d*u >= difmuk*difmuk*difmuk) return ignpoi; } ///////////////////////////////////////////////////////////////////// // STEP P. PREPARATIONS FOR STEPS Q AND H. // (RECALCULATIONS OF PARAMETERS IF NECESSARY) // .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7. // THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE // APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK. // C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION. /////////////////////////////////////////////////////////////////////// if (mu != muold) { muold = mu; omega = 0.3989423/s; b1 = 4.166667E-2/mu; b2 = 0.3*b1*b1; c3 = 0.1428571*b1*b2; c2 = b2-15.0*c3; c1 = b1-6.0*b2+45.0*c3; c0 = 1.0-b1+3.0*b2-15.0*c3; c = 0.1069/mu; } wflag = -1; if (g >= 0.0) { // 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN) kflag = 0; wflag = 0; } else { wflag = 1; } for (;;) { // STEP Q. QUOTIENT ACCEPTANCE (RARE CASE) if ( wflag > 0) { do { //////////////////////////////////////////////////////////////// // STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL // DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT' // (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.) //////////////////////////////////////////////////////////////// e = sexpo(); u = ranf(); u += (u-1.0); t = 1.8+fsign(e,u); } while (t <= -0.6744); ignpoi = (long) (mu+s*t); fk = (double)ignpoi; difmuk = mu-fk; // 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN) kflag = 1; } else { wflag = 1; } for (;;) { if (ignpoi < 10) { px = -mu; py = pow(mu,(double)ignpoi)/ *(fact+ignpoi); x = (0.5-difmuk)/s; xx = x*x; fx = -0.5*xx; fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0); if (kflag <= 0) { if(fy-u*fy <= py*exp(px-fx)) return ignpoi; break; } } else { //////////////////////////////////////////////////////////////// // CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION // A0-A7 FOR ACCURACY WHEN ADVISABLE // .8333333E-1=1./12. .3989423=(2*PI)**(-.5) //////////////////////////////////////////////////////////////// del = 8.333333E-2/fk; del -= (4.8*del*del*del); v = difmuk/fk; if (fabs(v) > 0.25) { px = fk*log(1.0+v)-difmuk-del; } else { px = fk*v*v*(((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0)-del; } py = 0.3989423/sqrt(fk); x = (0.5-difmuk)/s; xx = x*x; fx = -0.5*xx; fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0); if (kflag <= 0) { if(fy-u*fy <= py*exp(px-fx)) return ignpoi; break; } } // STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION) if (c*fabs(u) > py*exp(px+e)-fy*exp(fx+e)) {break;} else {return ignpoi;} } } }

long matvec::ignuin (  const long    low, const long    high )

Definition at line 1093 of file stat2.cpp. References ignuin(). Referenced by ignuin(), and matvec::DiscreteUniformDist::sample(). { #define maxnum 2147483647L if (low > high) throw exception(" ignu(low,high): low > hig"); // long random(); static long ignuin,ign,maxnow,range,ranp1; range = high-low; if (range > maxnum) throw exception(" ignuin(low,high): high too large"); if (low == high) { ignuin = low; return ignuin; } //////////////////////////////////////////////////////////////////// // Number to be generated should be in range 0..RANGE // Set MAXNOW so that the number of integers in 0..MAXNOW is an // integral multiple of the number in 0..RANGE //////////////////////////////////////////////////////////////////// ranp1 = range+1; maxnow = maxnum/ranp1*ranp1; S40: // ign = ignlgi()-1; ign = rand(); if(ign > maxnow) goto S40; ignuin = low+ign%ranp1; return ignuin; #undef maxnum }

void initial_BGarray (  BG *    v, unsigned    n, double    a )  [inline]

Definition at line 121 of file bg.h. Referenced by matvec::SparseBGMatrix::mv(), matvec::SparseBGMatrix::resize(), matvec::SparseBGMatrix::solve(), matvec::SparseBGMatrix::solvrow(), matvec::SparseBGMatrix::SparseBGMatrix(), and sweep(). { for (unsigned i = 0; i<n; i++) { v[i] = a; } }

int matvec::inva22 (  Matrix< double > &    a22 )  [static]

Definition at line 638 of file pop_mblup.cpp. References matvec::Session::epsilon, and SESSION. Referenced by matvec::Population::add_Gv_inv(), and matvec::Population::relv_inv(). { double tem,det; det = (a22[0][0]*a22[1][1] - a22[0][1]*a22[1][0]); if (fabs(det) < SESSION.epsilon) return 0; det = 1.0/det; tem=a22[0][0]; a22[0][0] = det*a22[1][1]; a22[1][1] = det*tem; a22[0][1] *= -det; a22[1][0] *= -det; return 1; }

doubleMatrix * matvec::KP (  Data *    D, doubleMatrix &    SKM, const std::string &    lpl, const std::string &    censor )

Definition at line 832 of file glmm2.cpp. References matvec::Data::datasheet, matvec::Data::field_index(), matvec::Data::num_rows(), matvec::Matrix< double >::resize(), matvec::Matrix< double >::sortby(), matvec::Matrix< double >::submat(), and warning(). { // 0=censored 1=uncensored; doubleMatrix ret,*retpt; doubleMatrix SKM; int ny; ny=D->num_rows(); doubleMatrix lpl_mat; lpl_mat.resize(ny,2); int lpl_col; if(-1 ==(lpl_col=D->field_index(lpl))) { warning("KP: Can't find lpl %s \n",lpl.c_str()); return(NULL); } for(int i=0;i<ny;i++){ lpl_mat[i][0]=D->datasheet[lpl_col][i].double_val(); } int status_col; if(-1 ==(status_col=D->field_index(censor))) { warning("KP: Can't find status %s \n",censor.c_str()); return(NULL); }; for(int i=0;i<ny;i++){ lpl_mat[i][1]=D->datasheet[status_col][i].double_val(); // std::cout << i << " " << lpl_mat[i][0] << " " << lpl_mat[i][1] << endl; } lpl_mat.sortby(Matrix<double>::COLUMN,0); // std::cout << "\n Input doubleMatrix col 1 = time col 2: 0=censored 1=uncensored " << lpl_mat; int j; float numdead; float numcensored; float lagdead; float lagcensored; doubleMatrix skm_dead; skm_dead.resize(ny,3); j=0; lagdead=0.; lagcensored=0.; for(int i=0;i<ny;i++){ if(i < ny-1 && lpl_mat[i][0] == lpl_mat[i+1][0]){ if(lpl_mat[i][1]==1){ lagdead++; } if(lpl_mat[i][1]==0){ lagcensored++; } } else{ if(lpl_mat[i][1]==1){ numdead=1.+lagdead; numcensored=lagcensored; lagdead=0.; lagcensored=0.; skm_dead[j][0]=lpl_mat[i][0]; skm_dead[j][1]=numdead; skm_dead[j][2]=numcensored; j++; } if(lpl_mat[i][1]==0.){ numdead=lagdead; numcensored=1.+lagcensored; lagdead=0.; lagcensored=0.; skm_dead[j][0]=lpl_mat[i][0]; skm_dead[j][1]=numdead; skm_dead[j][2]=numcensored; j++; } } } //std::cout << "\n SKM dead vector = " << skm_dead; skm_dead=skm_dead.submat(0,0,j,3); //std::cout << "\n SKM dead vector = " << skm_dead; //doubleMatrix SKM; SKM.resize(j+1,3); float tot_alive; tot_alive=ny; int k; k=0; for(int i=0;i<j+1;i++){ if(i == 0){ SKM[k][0]=0.; SKM[k][1]=1.; k++; } if(i > 0 && skm_dead[i-1][1] != 0.){ SKM[k][0]=skm_dead[i-1][0]; SKM[k][1]=SKM[k-1][1]*(1.-(skm_dead[i-1][1]/tot_alive)); SKM[k-1][2]=SKM[k][1]/SKM[k-1][1]; SKM[k-1][2]/=SKM[k][0]-SKM[k-1][0]; tot_alive=tot_alive-(skm_dead[i-1][1]+skm_dead[i-1][2]); k++; } if(i > 0 && skm_dead[i-1][1] == 0.){ tot_alive=tot_alive-(skm_dead[i-1][1]+skm_dead[i-1][2]); } } SKM[k-1][2]=SKM[k-2][2]; //for(int ii=0;ii<k;ii++) cout << SKM[ii][0] << " " << SKM[ii][1] << " " << SKM[ii][2] <<endl; SKM_ans=SKM.submat(0,0,k,3); return(&SKM_ans); }

unsigned lidx (  const unsigned    i, const unsigned    j, const unsigned    n )  [static]

Definition at line 33 of file model_blup.cpp. Referenced by matvec::Model::add_Ag_diag(), matvec::Model::add_Ig_diag(), matvec::Model::compute_rhs_diag_mme(), ginverse2(), and preconditioning(). { return ((n+n-j-1)*j+i+i)/2; }

void matvec::linmin (  double *    p, double *    xi, int    n, double &    fret )  [static]

Definition at line 124 of file min_nr_powell.cpp. References brent(), f1dim(), mnbrak(), ncom, pcom, TOL, and xicom. Referenced by nr_powell(). { int j; double xx,xmin,fx,fb,fa,bx,ax; ncom = n; if(n>0){ pcom = new double [n]; } else { pcom = 0; } if(n>0){ xicom = new double [n]; } else { xicom = 0; } for (j=0; j<n; j++) { pcom[j] = p[j]; xicom[j] = xi[j]; } ax=0.0; xx=1.0; bx=0.0; fa=0.0; fx=0.0; fb=0.0; xmin=0.0; mnbrak(ax,xx,bx,fa,fx,fb,f1dim); fret = brent(ax,xx,bx,f1dim,TOL,xmin); for (j=0; j<n; j++) { xi[j] *= xmin; p[j] += xi[j]; } if(xicom){ delete [] xicom; xicom=0; } if(pcom){ delete [] pcom; pcom=0; } }

template<class T> Matrix<T> log (  const Matrix< T > &    a )

Definition at line 1110 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::log); }

Field log (  Field &    a )

Definition at line 770 of file field.cpp. References matvec::Field::map(). {return a.map(std::log);}

BG log (  const BG &    a )

Definition at line 136 of file bg.cpp. References matvec::BG::D, matvec::BG::d, matvec::BG::f, and matvec::Matrix< double >::transpose(). Referenced by Beta_cdf(), betain(), matvec::GLMM::Build_SinMat(), ChiSquare_cdf(), ChiSquare_inv(), matvec::SparseBGMatrix::factor(), gamdev(), gammin(), gammln(), gasdev(), matvec::SparseBGMatrix::logdet(), Normal_inv(), sgamma(), t_cdf(), and t_expected_value(). { BG r; r.f = std::log(a.f); r.d = 1.0/a.f*a.d; r.D = 1.0/a.f*a.D - 1.0/(a.f*a.f) * a.d*a.d.transpose(); return r; }

template<class T> Matrix<T> log10 (  const Matrix< T > &    a )

Definition at line 1111 of file matrix.h. References matvec::Matrix< T >::apply(). { return a.apply(std::log10); }

Field log10 (  Field &    a )

Definition at line 771 of file field.cpp. References matvec::Field::map(). {return a.map(std::log10);}

void matvec::lubksb (  double **    a, int    n, int *    indx, double *    b )

Definition at line 389 of file ginverse1.cpp. References lubksb(). { int i,j; int ip, ii = -1; double sum; for (i=0;i<n;i++) { ip = indx[i]; sum = b[ip]; b[ip] = b[i]; if (ii>=0) { for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; } else if (sum) { ii=i; } b[i] = sum; } for (i=n-1; i>=0; i--) { sum = b[i]; for (j=i+1;j<n;j++) sum -= a[i][j]*b[j]; b[i] = sum/a[i][i]; } }

void matvec::lubksb (  BG **    a, int    n, int *    indx, BG *    b )

Definition at line 414 of file ginverse1.cpp. References lubksb(). Referenced by matvec::doubleMatrix::inv(), matvec::BGMatrix::inv(), matvec::doubleMatrix::lu_solve(), matvec::BGMatrix::lu_solve(), and lubksb(). { int i,j; int ip, ii = -1; BG sum; for (i=0;i<n;i++) { ip = indx[i]; sum = b[ip]; b[ip] = b[i]; if (ii>=0) { for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; } else if (sum!=0.0) { ii=i; } b[i] = sum; } for (i=n-1; i>=0; i--) { sum = b[i]; for (j=i+1;j<n;j++) sum -= a[i][j]*b[j]; b[i] = sum/a[i][i]; } }

void matvec::ludcmp (  double **    a, int    n, int *    indx, int &    d, double    tol )

Definition at line 291 of file ginverse1.cpp. References ludcmp(). { int i,j,k,imax; double big,dum,sum,temp; double *vv = new double [n]; d=1; for (i=0;i<n;i++) { big = 0.0; for (j=0;j<n;j++) if ((temp=std::fabs(a[i][j])) > big) big=temp; if (std::fabs(big) < tol) throw exception("ludcmp(): Singular matrix"); vv[i] = 1.0/big; } for (j=0;j<n;j++) { for (i=0;i<j;i++) { sum = a[i][j]; for (k=0;k<i;k++) sum -= a[i][k]*a[k][j]; a[i][j]=sum; } big=0.0; imax = j; for (i=j;i<n;i++) { sum = a[i][j]; for (k=0;k<j;k++) sum -= a[i][k]*a[k][j]; a[i][j] = sum; if ( (dum = vv[i]*std::fabs(sum)) >= big) { big = dum; imax = i; } } if (j != imax) { for (k = 0; k<n; k++) { dum = a[imax][k]; a[imax][k] = a[j][k]; a[j][k] = dum; } d = -d; vv[imax] = vv[j]; } indx[j] = imax; if (std::fabs(a[j][j]) < tol) throw exception("ludcmp(): Singular matrix"); dum = 1.0/(a[j][j]); for (i=j+1;i<n;i++) a[i][j] *= dum; } if(vv){ delete [] vv; vv=0; } }

void matvec::ludcmp (  BG **    a, int    n, int *    indx, int &    d, double    tol )

Definition at line 340 of file ginverse1.cpp. References ludcmp(). Referenced by matvec::doubleMatrix::det(), matvec::BGMatrix::det(), matvec::doubleMatrix::inv(), matvec::BGMatrix::inv(), matvec::doubleMatrix::lu_solve(), matvec::BGMatrix::lu_solve(), and ludcmp(). { int i,j,k,imax; BG big,dum,sum,temp; BG *vv = new BG [n]; d=1; for (i=0;i<n;i++) { big = 0.0; for (j=0;j<n;j++) if ((temp=fabs(a[i][j])) > big) big=temp; if (fabs(big) < tol) throw exception("ludcmp(): Singular matrix"); vv[i] = 1.0/big; } for (j=0;j<n;j++) { for (i=0;i<j;i++) { sum = a[i][j]; for (k=0;k<i;k++) sum -= a[i][k]*a[k][j]; a[i][j]=sum; } big=0.0; imax = j; for (i=j;i<n;i++) { sum = a[i][j]; for (k=0;k<j;k++) sum -= a[i][k]*a[k][j]; a[i][j] = sum; if ( (dum = vv[i]*fabs(sum)) >= big) { big = dum; imax = i; } } if (j != imax) { for (k = 0; k<n; k++) { dum = a[imax][k]; a[imax][k] = a[j][k]; a[j][k] = dum; } d = -d; vv[imax] = vv[j]; } indx[j] = imax; if (fabs(a[j][j]) < tol) throw exception("ludcmp(): Singular matrix"); dum = 1.0/(a[j][j]); for (i=j+1;i<n;i++) a[i][j] *= dum; } if(vv){ delete [] vv; vv=0; } }

void matvec::mnbrak (  double &    ax, double &    bx, double &    cx, double &    fa, double &    fb, double &    fc, double(*    func)(double) )  [static]

Definition at line 182 of file min_nr_powell.cpp. References GLIMIT, GOLD, MAX, SHFT, SIGN, and TINY. Referenced by linmin(). { double ulim,u,r,q,fu,dum; fa=(*func)(ax); fb=(*func)(bx); if (fb > fa) { SHFT(dum,ax,bx,dum) SHFT(dum,fb,fa,dum) } cx = bx+GOLD*(bx-ax); fc = (*func)(cx); while (fb > fc) { r=(bx-ax)*(fb-fc); q=(bx-cx)*(fb-fa); u=(bx)-((bx-cx)*q-(bx-ax)*r)/ (2.0*SIGN(MAX(fabs(q-r),TINY),q-r)); ulim = bx+GLIMIT*(cx-bx); if ((bx-u)*(u-cx) > 0.0) { fu=(*func)(u); if (fu < fc) { ax = bx; bx = u; fa = fb; fb = fu; return; } else if (fu > fb) { cx = u; fc = fu; return; } u= cx +GOLD*(cx-bx); fu=(*func)(u); } else if ((cx-u)*(u-ulim) > 0.0) { fu=(*func)(u); if (fu < fc) { SHFT(bx,cx,u,cx+GOLD*(cx-bx)) SHFT(fb,fc,fu,(*func)(u)) } } else if ((u-ulim)*(ulim-cx) >= 0.0) { u=ulim; fu=(*func)(u); } else { u=(cx)+GOLD*(cx-bx); fu=(*func)(u); } SHFT(ax,bx,cx,u) SHFT(fa,fb,fc,fu) } }

void nest_xaction (  std::string &    s )  [static]

Definition at line 2137 of file model.cpp. References nest_xaction(). Referenced by matvec::Model::build_model_term(), and nest_xaction(). { int i,j; int n = s.size(); for (i=0; i<n; i++) { if ( s[i] == '(' ) s[i] = '*'; if ( s[i] == ')' ) s[i] = ' '; } // replace space around * with * i = 0; while (i<n) { if (s[i] == '*' ) { j = i; while (s[--j] == ' ') s[j]='*'; while (s[++i] == ' ') s[i]='*'; } i++; } }

Chromosome* new_Chromosome_vec (  const unsigned    n )

Genome * matvec::new_Genome_vec (  const unsigned    m, GeneticDist *    D = 0 )

Definition at line 26 of file genome.cpp. References check_ptr(), and matvec::Genome::remodel(). Referenced by matvec::Population::build_pop_gamete(). { Genome *T = 0; if (m) { T = new Genome[m]; check_ptr(T); } if (D) { for (unsigned i=0; i<m; i++) T[i].remodel(D); } return T; }

HashNode* new_HashNode_vec (  const unsigned    n, const size_t    s )

RNormal* new_n_posterior_vec (  const int    m )

NuFamily* new_NuFamily_vec (  const unsigned    m )

long matvec::next_prime (  long    n )

Definition at line 103 of file hashtable.cpp. Referenced by matvec::HashTable::insert(), matvec::SparseMatrix::resize(), matvec::SparseBGMatrix::resize(), matvec::HashTable::resize(), matvec::SparseBGMatrix::SparseBGMatrix(), and matvec::SparseMatrix::SparseMatrix(). { long prime; // tentative prime long factor; // prime factor int is_prime; // 'prime' is a prime number n = (n < 0L) ? -n : n; // be safe -- force n positive prime = 2L*(n/2L) + 1L; // next odd number is_prime = (prime <= 3L); while (!is_prime) { factor = 5L; is_prime = (prime % 2L) && (prime % 3L); while (is_prime && (factor*factor <= prime)) { if ((prime % factor) == 0L) is_prime = 0; else if ((prime % (factor + 2L)) == 0L) is_prime = 0; else // factor+4 is divisible by 3 (every 3rd odd is) factor += 6L; } if (!is_prime) prime += 2L; } return (prime); }

int matvec::nonsymm_eigen (  const int    job, double **    A, const int    n, double *    rr, double *    ri, double **    vr, double **    vi, const int    maxiter = 50 )

Definition at line 590 of file nonsymm_eigen.cpp. References balance(), elemhess(), nonsymm_eigen_core(), sqrt(), and unbalance(). Referenced by matvec::doubleMatrix::eigen(). { ////////////////////////////////////////////////////////////////////////// // This nonsym_eigen() works for eigenvalue/vector analysis // for real general square matrix A // A will be destroyed // rr,ri are vectors containing eigenvalues // vr,vi are matrice containing eigenvectors // maxiter is max. no. of iterations to converge, usually 30 // note that A*v(:,j) = v(:,j)*r(j) where v(:,j) means jth column // job = 0, no eigenvectors are needed // 1, eigenvectors are also needed // // Algorithm: Handbook for Automatic Computation, vol 2 // by Wilkinson and Reinsch, 1971 // most of source codes were taken from a public domain // solftware called MATCALC. // Credits: goes to the authors of MATCALC // // return -2 out of memory // -1 not converged // 0 real eigenvalues/vectors // 1 complex eigenvalues/vectors // // Tianlin Wang at University of Illinois // Thu May 6 15:22:31 CDT 1993 ////////////////////////////////////////////////////////////////////////// double base = 2.0; // base of floating point arithmetic double digits = 53.0; // no. of digits to the base in the fraction double eps = pow(base,1.0-digits); double tiny = sqrt(eps); int low,hi,i,k=0; double *w; if(n>0){ w = new double [n]; } else { w = 0; } if (!w) return -2; balance(A,n,&low,&hi,w,base); elemhess(job,A,n,low,hi,vr,vi); if (-1 == nonsymm_eigen_core(job,A,n,low,hi,rr,ri,vr,vi,eps,maxiter)) { k = -1; } else { if (job) unbalance(n,vr,vi,low,hi,w); for (i=0; i<n; i++) if (fabs(ri[i]) > tiny) {k = 1; break;} } if(w){ delete [] w; w=0; } return k; }

int nonsymm_eigen_core (  const int    job, double **    mat, const int    n, const int    low, const int    hi, double *    valr, double *    vali, double **    vr, double **    vi, const double    eps, const int    maxiter )

Definition at line 262 of file nonsymm_eigen.cpp. References min, and sqrt(). Referenced by nonsymm_eigen(). { //////////////////////////////////////////////////////////////// // Calculate eigenvalues and eigenvectors of a real upper // Hessenberg matrix // Return 1 if converges successfully and 0 otherwise ///////////////////////////////////////////////////////////////// std::complex<double> v; double p,q,r,s,t,w,x,y,z,ra,sa,norm; int niter,en,i,j,k,l,m; for (i=0; i<n; i++) { valr[i]=0.0; vali[i]=0.0; } /* store isolated roots and calculate norm */ for (norm = 0,i = 0; i < n; i++) { for (j = std::max(0,i-1); j < n; j++) { norm += fabs(mat[i][j]); } if (i < low || i > hi) valr[i] = mat[i][i]; } t = 0; en = hi; while (en >= low) { niter = 0; for (;;) { /* look for single small subdiagonal element */ for (l = en; l > low; l--) { s = fabs(mat[l-1][l-1]) + fabs(mat[l][l]); if (s == 0) s = norm; if (fabs(mat[l][l-1]) <= eps * s) break; } /* form shift */ x = mat[en][en]; if (l == en) { /* one root found */ valr[en] = x + t; if (job) mat[en][en] = x + t; en--; break; } y = mat[en-1][en-1]; w = mat[en][en-1] * mat[en-1][en]; if (l == en - 1) { /* two roots found */ p = (y - x) / 2; q = p * p + w; z = sqrt(fabs(q)); x += t; if (job) { mat[en][en] = x; mat[en-1][en-1] = y + t; } if (q < 0) { /* complex pair */ valr[en-1] = x+p; vali[en-1] = z; valr[en] = x+p; vali[en] = -z; } else { /* real pair */ z = (p < 0) ? p - z : p + z; valr[en-1] = x + z; valr[en] = (z == 0) ? x + z : x - w / z; if (job) { x = mat[en][en-1]; s = fabs(x) + fabs(z); p = x / s; q = z / s; r = sqrt(p*p+q*q); p /= r; q /= r; for (j = en - 1; j < n; j++) { z = mat[en-1][j]; mat[en-1][j] = q * z + p * mat[en][j]; mat[en][j] = q * mat[en][j] - p*z; } for (i = 0; i <= en; i++) { z = mat[i][en-1]; mat[i][en-1] = q * z + p * mat[i][en]; mat[i][en] = q * mat[i][en] - p*z; } for (i = low; i <= hi; i++) { z = vr[i][en-1]; vr[i][en-1] = q*z + p*vr[i][en]; vr[i][en] = q*vr[i][en] - p*z; } } } en -= 2; break; } if (niter == maxiter) return(-1); if (niter != 0 && niter % 10 == 0) { t += x; for (i = low; i <= en; i++) mat[i][i] -= x; s = fabs(mat[en][en-1]) + fabs(mat[en-1][en-2]); x = y = 0.75 * s; w = -0.4375 * s * s; } niter++; /* look for two consecutive small subdiagonal elements */ for (m = en - 2; m >= l; m--) { z = mat[m][m]; r = x - z; s = y - z; p = (r * s - w) / mat[m+1][m] + mat[m][m+1]; q = mat[m+1][m+1] - z - r - s; r = mat[m+2][m+1]; s = fabs(p) + fabs(q) + fabs(r); p /= s; q /= s; r /= s; if (m == l || fabs(mat[m][m-1]) * (fabs(q)+fabs(r)) <= eps * (fabs(mat[m-1][m-1]) + fabs(z) + fabs(mat[m+1][m+1])) * fabs(p)) break; } for (i = m + 2; i <= en; i++) mat[i][i-2] = 0; for (i = m + 3; i <= en; i++) mat[i][i-3] = 0; /* double QR step involving rows l to en and columns m to en */ for (k = m; k < en; k++) { if (k != m) { p = mat[k][k-1]; q = mat[k+1][k-1]; r = (k == en - 1) ? 0 : mat[k+2][k-1]; if ((x = fabs(p) + fabs(q) + fabs(r)) == 0) continue; p /= x; q /= x; r /= x; } s = sqrt(p*p+q*q+r*r); if (p < 0) s = -s; if (k != m) { mat[k][k-1] = -s * x; } else if (l != m) { mat[k][k-1] = -mat[k][k-1]; } p += s; x = p / s; y = q / s; z = r / s; q /= p; r /= p; /* row modification */ for (j = k; j <= (!job ? en : n-1); j++){ p = mat[k][j] + q * mat[k+1][j]; if (k != en - 1) { p += r * mat[k+2][j]; mat[k+2][j] -= p * z; } mat[k+1][j] -= p * y; mat[k][j] -= p * x; } j = std::min(en,k+3); /* column modification */ for (i = (!job ? l : 0); i <= j; i++) { p = x * mat[i][k] + y * mat[i][k+1]; if (k != en - 1) { p += z * mat[i][k+2]; mat[i][k+2] -= p*r; } mat[i][k+1] -= p*q; mat[i][k] -= p; } if (job) { /* accumulate transformations */ for (i = low; i <= hi; i++) { p = x * vr[i][k] + y * vr[i][k+1]; if (k != en - 1) { p += z * vr[i][k+2]; vr[i][k+2] -= p*r; } vr[i][k+1] -= p*q; vr[i][k] -= p; } } } } } if (!job) return(0); if (norm != 0) { /* back substitute to find vectors of upper triangular form */ for (en = n-1; en >= 0; en--) { p = valr[en]; if ((q = vali[en]) < 0) { /* complex vector */ m = en - 1; if (fabs(mat[en][en-1]) > fabs(mat[en-1][en])) { mat[en-1][en-1] = q / mat[en][en-1]; mat[en-1][en] = (p - mat[en][en]) / mat[en][en-1]; } else { v = std::polar(0.0,-mat[en-1][en])/std::polar(mat[en-1][en-1]-p,q); mat[en-1][en-1] = v.real(); mat[en-1][en] = v.imag(); } mat[en][en-1] = 0; mat[en][en] = 1; for (i = en - 2; i >= 0; i--) { w = mat[i][i] - p; ra = 0; sa = mat[i][en]; for (j = m; j < en; j++) { ra += mat[i][j] * mat[j][en-1]; sa += mat[i][j] * mat[j][en]; } if (vali[i] < 0) { z = w; r = ra; s = sa; } else { m = i; if (vali[i] == 0) { v = std::polar(-ra,-sa)/std::polar(w,q); mat[i][en-1] = v.real(); mat[i][en] = v.imag(); } else { /* solve complex equations */ x = mat[i][i+1]; y = mat[i+1][i]; v = std::complex<double>((valr[i]- p)*(valr[i]-p) + vali[i]*vali[i] - q*q, (valr[i] - p)*2*q ); if ((fabs(v.real()) + fabs(v.imag())) == 0) { v = std::complex<double>(eps * norm * (fabs(w) + fabs(q) + fabs(x) + fabs(y) + fabs(z)), v.imag()); } v = std::polar(x*r-z*ra+q*sa,x*s-z*sa-q*ra) / v; mat[i][en-1] = v.real(); mat[i][en] = v.imag(); if (fabs(x) > fabs(z) + fabs(q)) { mat[i+1][en-1] = (-ra - w * mat[i][en-1] + q * mat[i][en]) / x; mat[i+1][en] = (-sa - w * mat[i][en] - q * mat[i][en-1]) / x; } else { v = std::polar(-r-y*mat[i][en-1], -s-y*mat[i][en]) / std::polar(z,q); mat[i+1][en-1] = v.real(); mat[i+1][en] = v.imag(); } } } } } else if (q == 0) { /* real vector */ m = en; mat[en][en] = 1; for (i = en - 1; i >= 0; i--) { w = mat[i][i] - p; r = mat[i][en]; for (j = m; j < en; j++) { r += mat[i][j] * mat[j][en]; } if (vali[i] < 0) { z = w; s = r; } else { m = i; if (vali[i] == 0) { if ((t = w) == 0) t = eps * norm; mat[i][en] = -r / t; } else { /* solve real equations */ x = mat[i][i+1]; y = mat[i+1][i]; q = (valr[i] - p) * (valr[i] - p) + vali[i]*vali[i]; t = (x * s - z * r) / q; mat[i][en] = t; if (fabs(x) <= fabs(z)) { mat[i+1][en] = (-s - y * t) / z; } else { mat[i+1][en] = (-r - w * t) / x; } } } } } } /* vectors of isolated roots */ for (i = 0; i < n; i++) { if (i < low || i > hi) { for (j = i; j < n; j++) { vr[i][j] = mat[i][j]; } } } /* multiply by transformation matrix */ for (j = n-1; j >= low; j--) { m = std::min(j,hi); for (i = low; i <= hi; i++) { for (z = 0,k = low; k <= m; k++) { z += vr[i][k] * mat[k][j]; } vr[i][j] = z; } } } /* rearrange complex eigenvectors */ for (j = 0; j < n; j++) { if (vali[j] != 0) { for (i = 0; i < n; i++) { vi[i][j] = vr[i][j+1]; vr[i][j+1] = vr[i][j]; vi[i][j+1] = -vi[i][j]; } j++; } } return(0); }

double matvec::Normal_cdf (  const double    x )

Definition at line 163 of file stat1.cpp. References erf(), and sqrt(). Referenced by matvec::NormalDist::cdf(), gammin(), t_cdf(), and matvec::LogNormalDist::variance(). { // ****************************************************************** // return the cumulative distribution value for a standard normal // ***************************************************************** double pr = 0.5; if (x == 0.0) return pr; register double c = 1.0/sqrt(2.0); // c = 1/sqrt(2) if (x <0.0) { pr = 0.5 - 0.5*erf(-x*c);} else { pr = 0.5 + 0.5*erf(x*c);} return pr; }

double matvec::Normal_inv (  const double    p )

Definition at line 571 of file stat1.cpp. References log(), and sqrt(). Referenced by ChiSquare_inv(), matvec::LogNormalDist::inv(), and matvec::NormalDist::inv(). { ////////////////////////////////////////////////////////////////// // ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3 // // Produces the standard normal deviate Z corresponding to a given // lower tail area of P; Z is accurate to about 1 part in 10**16. ////////////////////////////////////////////////////////////////// if (p <= 0.0 || p >= 1.0 ) throw exception(" Normal_inv(): bad arg, value in (0,1) expected"); if (p == 0.5) return 0.0; double *a = new double [8]; double *b = new double [8]; double *c = new double [8]; double *d = new double [8]; double *e = new double [8]; double *f = new double [8]; a[0] = 3.3871328727963666080e+0; a[1] = 1.3314166789178437745e+2; a[2] = 1.9715909503065514427e+3; a[3] = 1.3731693765509461125e+4; a[4] = 4.5921953931549871457e+4; a[5] = 6.7265770927008700853e+4; a[6] = 3.3430575583588128105e+4; a[7] = 2.5090809287301226727e+3; b[0] = 0.0; b[1] = 4.2313330701600911252e+1; b[2] = 6.8718700749205790830e+2; b[3] = 5.3941960214247511077e+3; b[4] = 2.1213794301586595867e+4; b[5] = 3.9307895800092710610e+4; b[6] = 2.8729085735721942674e+4; b[7] = 5.2264952788528545610e+3; c[0] = 1.42343711074968357734e+0; c[1] = 4.63033784615654529590e+0; c[2] = 5.76949722146069140550e+0; c[3] = 3.64784832476320460504e+0; c[4] = 1.27045825245236838258e+0; c[5] = 2.41780725177450611770e-1; c[6] = 2.27238449892691845833e-2; c[7] = 7.74545014278341407640e-4; d[0] = 0.0; d[1] = 2.05319162663775882187e+0; d[2] = 1.67638483018380384940e+0; d[3] = 6.89767334985100004550e-1; d[4] = 1.48103976427480074590e-1; d[5] = 1.51986665636164571966e-2; d[6] = 5.47593808499534494600e-4; d[7] = 1.05075007164441684324e-9; e[0] = 6.65790464350110377720e+0; e[1] = 5.46378491116411436990e+0; e[2] = 1.78482653991729133580e+0; e[3] = 2.96560571828504891230e-1; e[4] = 2.65321895265761230930e-2; e[5] = 1.24266094738807843860e-3; e[6] = 2.71155556874348757815e-5; e[7] = 2.01033439929228813265e-7; f[0] = 0.0; f[1] = 5.99832206555887937690e-1; f[2] = 1.36929880922735805310e-1; f[3] = 1.48753612908506148525e-2; f[4] = 7.86869131145613259100e-4; f[5] = 1.84631831751005468180e-5; f[6] = 1.42151175831644588870e-7; f[7] = 2.04426310338993978564e-15; double split1 = 0.425; double split2 = 5.0; double const1 = 0.180625; double const2 = 1.6; double q,r,w1,w2,ppnd16; q = p - 0.5; if (fabs(q) <= split1) { r = const1 - q*q; w1 = ((((((a[7]*r+a[6])*r+a[5])*r+a[4])*r+a[3])*r+a[2])*r+a[1])*r+a[0]; w2 = ((((((b[7]*r+b[6])*r+b[5])*r+b[4])*r+b[3])*r+b[2])*r+b[1])*r+1.0; ppnd16 = q*w1/w2; } else { if (q < 0.0) {r = p;} else {r = 1.0-p;} if (r <= 0.0) { if(a){ delete [] a; a=0; } if(b){ delete [] b; b=0; } if(c){ delete [] c; c=0; } if(d){ delete [] d; d=0; } if(e){ delete [] e; e=0; } if(f){ delete [] f; f=0; } return 0.0; } r = sqrt(-log(r)); if (r <= split2) { r -= const2; w1= ((((((c[7]*r+c[6])*r+c[5])*r+c[4])*r+c[3])*r+c[2])*r+c[1])*r+c[0]; w2= ((((((d[7]*r+d[6])*r+d[5])*r+d[4])*r+d[3])*r+d[2])*r+d[1])*r+1.0; } else { r -= split2; w1= ((((((e[7]*r+e[6])*r+e[5])*r+e[4])*r+e[3])*r+e[2])*r+e[1])*r+e[0]; w2= ((((((f[7]*r+f[6])*r+f[5])*r+f[4])*r+f[3])*r+f[2])*r+f[1])*r+1.0; } ppnd16 = w1/w2; if (q< 0.0) ppnd16 = -ppnd16; } if(a){ delete [] a; a=0; } if(b){ delete [] b; b=0; } if(c){ delete [] c; c=0; } if(d){ delete [] d; d=0; } if(e){ delete [] e; e=0; } if(f){ delete [] f; f=0; } return ppnd16; }

void matvec::normal_loginClasses (  matvec::Observe &    observe )

Definition at line 1749 of file glmm1.cpp. References matvec::Observe::estimate, matvec::doubleMatrix::ginv1(), matvec::Observe::H, matvec::Observe::like_val, matvec::doubleMatrix::logdet(), matvec::doubleMatrix::mat_exp(), matvec::doubleMatrix::mat_exp_der(), matvec::doubleMatrix::mat_log(), matvec::Observe::numtrait, matvec::Observe::nvc, matvec::Observe::pattern, matvec::doubleMatrix::quadratic(), matvec::Observe::rec, matvec::Observe::Resid_Sin_Mat, matvec::Observe::Resid_Var, matvec::Observe::resid_var, matvec::Observe::residual, matvec::Vector< T >::resize(), matvec::Matrix< double >::resize(), matvec::Observe::rpy, matvec::Observe::rve, matvec::Matrix< double >::transpose(), matvec::Observe::varcomp, and matvec::Observe::y. Referenced by matvec::GLMM::normalLog(). { int numtrait=observe.numtrait; doubleMatrix H,variance; doubleMatrix Alog; doubleMatrix P,R; Vector<double> D; // double final=0; //if(observe.param) final=*((double *) observe.param[0]); int k,t1,t2; if(observe.rec < 0) { if(observe.nvc) { if(observe.nvc != (numtrait*(numtrait+1))/2) { throw matvec::exception("Error: Should be normal(,%d)"); return; } //cout << variance << observe.resid_var[0]; *observe.resid_var=observe.resid_var[0]; Alog=(observe.resid_var[0]).mat_log(); for(k=0,t1=0;t1<numtrait;t1++) for(t2=t1;t2<numtrait;t2++,k++) observe.varcomp[0][k]=Alog[t1][t2]; } return; } Vector<double> mean,y; Alog.resize(numtrait,numtrait); Vector<doubleMatrix> part(observe.nvc); for(k=0;k<observe.nvc;k++) part[k].resize(numtrait,numtrait); if(observe.nvc) { for(k=0,t1=0;t1<numtrait;t1++) for(t2=t1;t2<numtrait;t2++,k++) { Alog[t1][t2]=observe.varcomp[0][k]; Alog[t2][t1]=observe.varcomp[0][k]; } variance=Alog.mat_exp(); //cout << variance; Alog=variance.mat_log(); for(k=0,t1=0;t1<numtrait;t1++) for(t2=t1;t2<numtrait;t2++,k++) observe.varcomp[0][k]=Alog[t1][t2]; *observe.resid_var=variance; variance.mat_exp_der(part); doubleMatrix Atmp; /* for(int ii=0,t1=0;t1 <numtrait;t1++) for(t2=t1;t2<numtrait;t2++,ii++) { Atmp=Alog; Atmp[t1][t2]+=.0001; Atmp[t2][t1]=Atmp[t1][t2]; Atmp=(Atmp.mat_exp()-Alog.mat_exp())/.0001; cout << ii << part[ii]<< endl; cout << ii << Atmp << endl; } */ } else{ variance=observe.resid_var[0]; } //cout << "\n Variance " << variance << *observe.resid_var; doubleMatrix Miss; Miss.resize(numtrait,numtrait); //cout << "\n Pattern" ; for(t1=0;t1<numtrait;t1++) { //cout <<' ' <<observe.pattern[t1]; if(observe.pattern[t1]=='1') Miss[t1][t1]=1; } //Miss.identity(); //static double mean,H,variance,scale,y; y=observe.y; y=Miss*y; //cout << "\n y\n" << y; // if(observe.estimate[0] < -5.) observe.estimate[0]=-5.; // if(observe.estimate[0] > 5.) observe.estimate[0]=5.; mean=Miss*observe.estimate; variance=*observe.resid_var; variance=Miss*variance*Miss; //cout << "\n Variance \n" << variance; observe.Resid_Var=variance; variance=variance.ginv1(); //cout << "\n Variance inv \n" << variance; observe.like_val=.5*(variance.logdet()-variance.quadratic((y-mean),(y-mean))) ; observe.H=Miss; H=observe.H; //observe.Hinv[0][0]=1./H; observe.rve=H.transpose()*variance*H; observe.rpy=H.transpose()*variance*y; //cout <<"\n rve rpy "<<observe.rve << " " << observe.rpy<< " H "<< H << "pattern" << observe.pattern <<" y"<< y <<" var" << variance; observe.residual=y-mean; if(observe.nvc) { for(k=0,t1=0;t1<numtrait;t1++) for(t2=t1;t2<numtrait;t2++,k++) { observe.Resid_Sin_Mat[k]=variance*Miss*part[k]*Miss*variance; // cout << k << " " << part[k] << observe.Resid_Sin_Mat[k] << endl; } } }

double matvec::nr_powell (  Data *    D, Model *    M, Vector< double > &    p, Matrix< double > &    xi, int    nvc, int &    iter, double    ftol )

Definition at line 63 of file min_nr_powell.cpp. References matvec::Vector< T >::begin(), funk(), linmin(), nfunk, SQR, and warning(). Referenced by matvec::Model::vce_dfreml(). { myD = D; myM = M; nfunk = 0; int i,ibig,j; double t,fptt,fp,del,fret; int maxiter = iter; Vector<double> pt(n); Vector<double> ptt(n); Vector<double> xit(n); fret = funk(p); for (j=0; j<n; j++) pt[j]=p[j]; for (iter=1; ;iter++) { fp = fret; ibig = 0; del = 0.0; for (i=0; i<n; i++) { for (j=0; j<n; j++) xit[j] = xi[j][i]; fptt = fret; linmin(p.begin(),xit.begin(),n,fret); if (fabs(fptt-fret) > del) { del = fabs(fptt-fret); ibig = i; } } if (2.0*fabs(fp-fret) <= ftol*(fabs(fp)+fabs(fret))) { iter = nfunk; return fret; } if (iter == maxiter) { warning(" Too many iterations in routine POWELL\n" "%d out of %d",iter,maxiter); return fret; } for (j=0; j<n; j++) { ptt[j] = 2.0*p[j] - pt[j]; xit[j] = p[j] - pt[j]; pt[j] = p[j]; } fptt = funk(ptt); if (fptt < fp) { t = 2.0*(fp-2.0*fret+fptt)*SQR(fp-fret-del)-del*SQR(fp-fptt); if (t < 0.0) { linmin(p.begin(),xit.begin(),n,fret); for (j=0; j<n; j++) xi[j][ibig] = xit[j]; } } } }

SparseMatrix operator * (  const double    s, const SparseMatrix &    A )

Definition at line 499 of file sparsematrix.cpp. References operator *(). {return operator*(A,s);}

SparseMatrix operator * (  const SparseMatrix &    A, const double    s )

Multiply SparseMatrix object by s. Definition at line 475 of file sparsematrix.cpp. References matvec::SparseMatrix::aav, matvec::SparseMatrix::dim, matvec::SparseMatrix::hash_srt, matvec::SparseMatrix::hsize, matvec::SparseMatrix::iiv, matvec::SparseMatrix::jjv, matvec::SparseMatrix::max_nz, matvec::SparseMatrix::nonzero, and matvec::SparseMatrix::srt_hash. { if ( A.dim == 0 ) throw exception("SparseMatrix is empty"); unsigned new_dim = A.dim; unsigned new_nz = A.nonzero; SparseMatrix B(new_dim,A.max_nz); B.nonzero = new_nz; memcpy(B.iiv,A.iiv,sizeof(int)*A.hsize); memcpy(B.jjv,A.jjv,sizeof(int)*A.hsize); if(A.hsize){ memcpy(B.hash_srt,A.hash_srt,sizeof(int)*(A.hsize+1)); memcpy(B.srt_hash,A.srt_hash,sizeof(int)*(A.hsize+1)); } double *A_pt = A.aav; double *New_pt = B.aav; for (unsigned i=0; i<A.hsize; i++) { if (A.iiv[i] > 0) { New_pt[i] = A_pt[i]*s; } } return B; }

Vector<double> operator * (  const SparseMatrix &    A, const Vector< double > &    v )

SparseMatrix times vector operation. Definition at line 395 of file sparsematrix.cpp. References matvec::SparseMatrix::aav, matvec::Vector< T >::begin(), matvec::SparseMatrix::dim, matvec::SparseMatrix::iiv, matvec::SparseMatrix::jjv, matvec::SparseMatrix::nonzero, matvec::Vector< T >::reserve(), matvec::Vector< T >::size(), and matvec::SparseMatrix::srt_hash. { if ( A.dim == 0) throw exception("SparseMatrix is empty"); if (A.dim != v.size()) throw exception("SparseMatrix*Vector<double>, size not conformable"); Vector<double> vec_out; int n = A.dim; int *A_ia = A.iiv-1; /* offset by 1 */ int *A_ja = A.jjv-1; double *A_a = A.aav-1; double *v_ve = v.begin()-1; double tmp,x; unsigned k,i,j; vec_out.reserve(n); double *out = vec_out.begin(); out--; // offset by 1 int ksrt=0; //SDK for(ksrt=0;ksrt<A.nonzero;ksrt++) { //SDK k=A.srt_hash[ksrt]; //SDK // for (k=1; k<=A.hsize; k++) { //SDK if (i = A_ia[k]) { j = A_ja[k]; x = A_a[k]; out[i] = out[i] + x*v_ve[j]; if (j>i) { out[j] = out[j] + x*v_ve[i]; } } } return vec_out; }

SparseBGMatrix operator * (  const BG    s, const SparseBGMatrix &    A )

Definition at line 503 of file sparsebgmatrix.cpp. References operator *(). {return operator*(A,s);}

SparseBGMatrix operator * (  const SparseBGMatrix &    A, const BG    s )

Multiply SparseBGMatrix object by s. Definition at line 479 of file sparsebgmatrix.cpp. References matvec::SparseBGMatrix::aav, matvec::SparseBGMatrix::dim, matvec::SparseBGMatrix::hash_srt, matvec::SparseBGMatrix::hsize, matvec::SparseBGMatrix::iiv, matvec::SparseBGMatrix::jjv, matvec::SparseBGMatrix::max_nz, matvec::SparseBGMatrix::nonzero, and matvec::SparseBGMatrix::srt_hash. { if ( A.dim == 0 ) throw exception("SparseBGMatrix is empty"); unsigned new_dim = A.dim; unsigned new_nz = A.nonzero; SparseBGMatrix B(new_dim,A.max_nz); B.nonzero = new_nz; memcpy(B.iiv,A.iiv,sizeof(int)*A.hsize); memcpy(B.jjv,A.jjv,sizeof(int)*A.hsize); if(A.hsize){ memcpy(B.hash_srt,A.hash_srt,sizeof(int)*(A.hsize+1)); memcpy(B.srt_hash,A.srt_hash,sizeof(int)*(A.hsize+1)); } BG *A_pt = A.aav; BG *New_pt = B.aav; for (unsigned i=0; i<A.hsize; i++) { if (A.iiv[i] > 0) { New_pt[i] = A_pt[i]*s; } } return B; }

Vector<BG> operator * (  const SparseBGMatrix &    A, const Vector< BG > &    v )

SparseBGMatrix times vector operation. Definition at line 398 of file sparsebgmatrix.cpp. References matvec::SparseBGMatrix::aav, matvec::Vector< T >::begin(), matvec::SparseBGMatrix::dim, matvec::SparseBGMatrix::iiv, matvec::SparseBGMatrix::jjv, matvec::SparseBGMatrix::nonzero, matvec::Vector< T >::reserve(), matvec::Vector< T >::size(), and matvec::SparseBGMatrix::srt_hash. { if ( A.dim == 0) throw exception("SparseBGMatrix is empty"); if (A.dim != v.size()) throw exception("SparseBGMatrix*Vector<BG>, size not conformable"); Vector<BG> vec_out; int n = A.dim; int *A_ia = A.iiv-1; /* offset by 1 */ int *A_ja = A.jjv-1; BG *A_a = A.aav-1; BG *v_ve = v.begin()-1; BG tmp,x; unsigned k,i,j; vec_out.reserve(n); BG *out = vec_out.begin(); out--; // offset by 1 int ksrt=0; //SDK for(ksrt=0;ksrt<A.nonzero;ksrt++) { //SDK k=A.srt_hash[ksrt]; //SDK // for (k=1; k<=A.hsize; k++) { //SDK if (i = A_ia[k]) { j = A_ja[k]; x = A_a[k]; out[i] = out[i] + x*v_ve[j]; if (j>i) { out[j] = out[j] + x*v_ve[i]; } } } return vec_out; }

template<class T> Matrix<T> operator * (  const T &    a, const Matrix< T > &    b )

Definition at line 1123 of file matrix.h. { return b * a; }

template<class T> Vector<T> operator * (  const Vector< T > &    a, const Matrix< T > &    b )  [inline]

Definition at line 1012 of file matrix.h. References matvec::Matrix< T >::begin(), matvec::Matrix< T >::num_cols(), matvec::Matrix< T >::num_rows(), and matvec::Vector< T >::size(). { if (a.size() != b.num_rows()) throw exception("Matrix<T>::operator*: size not conformable"); Vector<T> temp(b.num_cols()); T** me = b.begin(); T x; for (int i,j=0; j<b.num_cols(); ++j) { x = T(); for (i=0; i<b.num_rows(); ++i) x += a[i] * me[i][j]; temp[j] = x; } return temp; }

doubleMatrix matvec::operator * (  Vector< double > &    a, Vector< double > &    b )

Definition at line 1084 of file glmm2.cpp. References matvec::Vector< T >::size(). { int row,col; row=a.size(); col=b.size(); doubleMatrix temp(row,col); for(int i=0;i<row;i++) { for(int j=0;j<col;j++) { temp[i][j]=a[i]*b[j]; } } return temp; }

Field operator * (  const double    s, const Field &    B )

Definition at line 632 of file field.cpp. References matvec::Field::col_struct, matvec::Field::dat_vec, matvec::DataNode::double_val(), matvec::DataNode::missing, matvec::Field::ne, matvec::FieldStruct::type(), and warning(). { if (B.col_struct.type()=='S') { warning("s * Field: failed for string column"); return Field(); } DataNode *temp; if(B.ne>0){ temp = new DataNode [B.ne]; } else { temp = 0; } for (unsigned i=0; i<B.ne; i++) { if (B.dat_vec[i].missing) { temp[i].missing = 1; } else { temp[i].double_val( s * B.dat_vec[i].double_val()); } } return Field(B.ne,temp,B.col_struct); }

Field operator * (  const Field &    A, const double    s )

Definition at line 536 of file field.cpp. References matvec::Field::col_struct, matvec::Field::dat_vec, matvec::DataNode::double_val(), matvec::DataNode::missing, matvec::Field::ne, matvec::FieldStruct::type(), and warning(). { if (A.col_struct.type()=='S') { warning("Field * s: failed for string column"); return Field(); } DataNode *temp; if(A.ne>0){ temp = new DataNode [A.ne]; } else { temp = 0; } for (unsigned i=0; i<A.ne; i++) { if (A.dat_vec[i].missing) { temp[i].missing = 1; } else { temp[i].double_val(A.dat_vec[i].double_val() * s); } } return Field(A.ne,temp,A.col_struct); }

Field operator * (  const Field &    A, const Field &    B )

Definition at line 428 of file field.cpp. References matvec::Field::col_struct, matvec::Field::dat_vec, matvec::DataNode::double_val(), matvec::DataNode::missing, matvec::Field::ne, matvec::FieldStruct::type(), and warning(). { unsigned n = A.ne; if (n != B.ne) { warning("Col * Col: size not conformable"); return Field(); } if (A.col_struct.type()=='S' || B.col_struct.type()=='S') { warning("Col * Field: failed for string column"); return Field(); } DataNode *temp; if(n>0){ temp = new DataNode [n]; } else { temp = 0; } for (unsigned i=0; i<n; i++) { if (A.dat_vec[i].missing || B.dat_vec[i].missing) { temp[i].missing = 1; } else { temp[i].double_val(A.dat_vec[i].double_val()* B.dat_vec[i].double_val()); } } return Field(n,temp,A.col_struct); }

BG operator * (  const BG &    a, double    b )

Definition at line 189 of file bg.cpp. { BG r = b * a; return r; }

BG operator * (  double    a, const BG &    b )

Definition at line 168 of file bg.cpp. References matvec::BG::D, matvec::BG::d, and matvec::BG::f. { BG r; r.f = a*b.f; r.d = a*b.d; r.D = a*b.D; return r; }

BG operator * (  const BG &    a, const BG &    b )

Definition at line 99 of file bg.cpp. References matvec::BG::D, matvec::BG::d, matvec::BG::f, and matvec::Matrix< double >::transpose(). Referenced by operator *(). { BG r; r.f = a.f * b.f; r.d = a.f*b.d + a.d*b.f; r.D = a.f*b.D + a.d*b.d.transpose() + b.d*a.d.transpose() + a.D*b.f; return r; }

BG& operator *= (  BG &    a, double    b )

Definition at line 229 of file bg.cpp. { a = a * b; return a; }

<