*************************************************************************
*
* Changes in BLUPF90 and related programs
*
*************************************************************************
-----------------------------------------------------------------
Changes from Nov 2003 are inversely chronological (newest on top)
-----------------------------------------------------------------
Sept 25, 2007
Too long fields that used readln, nums and nums2 generated error messages
because bounds were not checked.
Increased input line length from 200 to 400 charaters.
Corrected statistics in renumf90 on covariables.
Done by I. Aguillar
July 3, 2007
Change of May 24 reverted due to slow hashing with sparse matrices.
May 25, 2007
Two additional operations added to FSPAK and FSPAK90: mutiplication by a
factor and solutions with the factor as LHS; programming by Juan Pablo Sanchez.
May 24, 2007
Storing elements in sparse_ija formt when the matrix was relatively dense
resulted in long computing time and large memory requirements. Prime
numbers in hashvr() were changed for bigger numbers. Modifications per
Andres Legarra.
April 9, 2007
In search.f90, sorting large alpha fields did not work with some
compilers, and positioning for the last element returned the element
next to the last.
February 26, 2007
Sorting in search.f90 returned incorrect results with multiple
sorting fields.
Feb 2, 2007
fspak-metis modified to remove an error message with a compiler.
Aug 1, 2006
Renumf90 generated cryptic messages when UPG were specified in_pedigrees and
PED_DEPTH was not zero.
July 3, 2006
Gibbs1 and later programs returned incorrect results under IBM AIX. For a
quick fix, rename file gibbsibm.f90 as gibbs.f90 in libs.
May 11, 2006
RENUMF90 stops with an error message if UPG is specified and the position
of year of birth is missing
Feb 2, 2006
CBLUP90IOD modified for improved stability; all thresholds are now
estimated; previously last threshold was set to 0.
Dec 10, 2005.
RENUMF90 now accepts a COMBINE keyword that allows a combination of
fields to be treated as a new field.
Sept 8, 2005
A new version of fspak.f is available to work with METIS, a rapid
ordering software. For details, read Readme.fspak-metis in the libs
directory.
Feb 8, 2005
Subroutine iounf is now overloaded to read/write integer and real*8 vectors.
Feb 7, 2004
Renumf90 modified for PED_DEPTH 3 and flexible output format allowing
for > 9,999,999 levels; also memory leak eliminated.
blup90iod modified to correctly print levels greater than 9,999,999.
Jan 5, 2005
Added version of Postgibbsf90 for SUSE64; does not seem to work under
SUSE32.
December 21, 2004
Renumf90 ocassionally produced an error message "pedigrees repeated".
October 28, 2004
Solutions for covariables missing for some traits were not zero
with some compilers. Subroutine find_addresses were replaced in
all programs except MRF90.
May 20, 2004
Presence of weight variable caused residual variance to
approach zero in all gibbs sampling programs.
April 13
Overloaded assignments (e.g., xx_ija=xx) caused memory leak
and subsequently "out of memory" errors with some
compilers including Lahey. This was corrected by changing
intent(out) to intent(inout) in subroutines implementing
those assignments.
February 20, 2004
Memory leak in gibbs2f90 and gibbs3f90 eliminated.
Weights in airemlf90 now function correctly.
Jan 24, 2004
Blup90iod was modified to implement Lidauer et al acceleration. Speed-up
depends on p = number of effects x number of traits; speedup
is small for small p but could be 5 times for p=100.
Dec 3
Conversion from hash to ija fromats caused deallocation error
messages mainly on HP and IBM/AIX workstations. Nullifying a
pointer in sparse.f90 fixed the problem.
Dec 2
Added to SPARSEM overloaded routines for copying matrix structures
of the same kind.
Overloaded subroutine RESET now checks whether pointers are
associated before deallocation. If they are not associated, there is
no deallocation.
Nov 14
AIREMLRES added to the distribution
-----------------------------------------------------------------
Changes below are in chronological order (newest on bottom)
-----------------------------------------------------------------
May 12/99
GIBBSF90 worked incorrectly when some variances were not estimated and
were set to 0. Fixed random number generators in rndgen.f90 so that
variances of 0 do not return NaN.
April-May 1999
Added matrix-scalar and matrix-vector multiplication to SPARSEM.
Added printing of rectangular matrices in DENSEOP.
Added a function diag that changes matrix to vector and vice versa in
DENSEOP.
Added subroutine pos_def that checks whether a matrix is semipositive-
definite, and if not, corrects it; in DENSEOP.
Added checking for semi-positive dispersion matrices in BLUPS1.F90;
also additional checks in REMLF90.
June-July, 1999
Fixed denseop.f90 so that it would work with the Fujitsu compiler on
Sun workstations.
Added AIREMLF90 by Shogo Tsuruta that makes REMLF90 converge 10-20 times
faster; the current version does not work yet with correlated effects.
Refined GIBBSF90 so that it seems to work with all problems except
very small examples.
September, 1999
Fixed a bug in GIBBSF90 where covariances set to zero in
initial covariances matrices were still estimated.
Module PCG added that computes solutions by Preconditioned Conjugate Gradient
with diagonal and non-diagonal preconditioners.
Made PCG the deafult solver in BLUPF90. . The convergence
is now faster, by a factor of 3-50 depending on model.
The whole package no longer works with real*4 precision, when rh=r4
in kind.f90. Fixes are planned.
October 1999
GIBBSF90 no longer updates covariances that are 0 originally. Same
correction needs to be applied to REMLF90i when covariances are 0
but all corresponding variances are not zero.
November 1999
Fixed crashes with IOUNF ocuuring with certain compilers where
preopened units had to be closed before they could be opened by IOUNF.
November 22, 1999
(log)determinant in FSPAK is now computed correctly with matrices not of
full rank.
November 26, 1999
In DENSEM, the Cholesky subroutines now create lower diagonal matrices
rather than symmetric matrices. Subroutine rndgen.f in directory GIBBS
modified to account for this.
January 26, 2000
Fixed crashes with addm() to large SPARSE_HASHM matrices
February 1, 2000
On HP, pointer arrays may not be nullified on creation by the system, causing
crashes. With two identical models except that the animal effect
was first in the first model and was last in the second model, remlf90
work fine with the first one and crashed with the second one. Partially fixed
by inserting
init(ainv(i))
before
call zero(ainv(i),n)
in one but not all places. Nullifying pointers on creation is a de-facto but
not the real standard.
February 8, 2000
Subroutine gen_invwishart destroyed the variance parameter.
February 24, 2000
Made subroutines reading and writing the parameter file part of module
textop. Now, all programs using these subroutines need to have
use textop
next to other "use" statements.
The module includes a new procedure nums2 that tokenizes a line into
numbers and strings separated by spaces. As opposed to the previous version
nums(), it can read real numbers.
Also included is a subroutine that can read optional Parameters from the
parameter file. The optional parameters are at the end of the parameter
file and have the form:
OPTION name str1 str2 ...
Description of these subroutines can be found as comments to
those subroutines in libs/blups1.f90
March 21
The likelihood in REMLF90 is now computed correctly and takes into account
not only G0 but also A so that comparisons of different models are possible.
It is not clear whether the likelihood constants are correct in models
with unknown parent groups and models with different random effects for
each trait.
April 13
Likelihood in REMLF90 was incorrect with some combination of missing trait.
Variance components were incorrect when weights were used.
Nested regression did not work if their number was < number of effects or
larger than any effect or trait column.
April 24
Took the gibbs.f90 module out of directory gibbs and moved to directory libs.
Unit prob.f90 used for calculation of probabilities/cdf in threshold models
was made a module and also moved to directory libs. The two programs were
arranged then logically, where all statistical routines were moved to prob.f90
and all gibss sampling manipulation routines were retained in gibbs.f90. Both
programs are now compiled into library stat.a. Appropriate changes were made
to makefiles and sources programs in directories gibbs and thr*. The
probability routines were updated with some for truncated normal distributions
useful for implementation of gibbs sampling in threshold univariate and
bivariate models.
May 2
Files prob.f90 and rndgen.f are now combined as prob.f90. Uniform random number
generator has been improved. It is also overloaded to operate for integer
arguments and varying intervals. Generator for truncated distribution has been
generalized. Argument "seed" has been removed from all functions. Seed can be
changed using a new function. The manual for module stat.a is Readme in
directory gibbs.
May 18
Residuals in REMLF90 became zero when trait values were very large. Possibly,
message about R being not positive-definite was printed. One of the fixes
involved changing the operational zero in module denseop from 10d-6
to 10d-10. If residuals in R are estimated as 0 in round 2, and observations
have very large values, scale the observations so that they are smaller.
June 12
Compilation of programs under Service Pack 4 of Absoft Linux compiler version
6.1 (or 6.2?) and older SGI generated an error message due to
optimization bugs. This has been corrected by a slight change in the
store_solutions subroutine.
June 14
Under sire model, the hash matrix would enlarge when it is less than 10%
full, resulting in very large memory requirements.
September 1
quad_form for SPARSE_HASH format generated accesses to 0 indices, which
caused crashes under DEC computers.
In prob.f90, gen_uniform did not sample randomly from all 100 elements,
and sometimes indices of 0 were generated causing random crashes.
In fspak90.f90, second routine has been added an external attribute,
which is required by some compilers.
In Gibbs1f90, df_random was not zeroed which caused strange message under
DEC for some data files.
September 13
New Gibbs sampling program available in gibbsfast; It creates only
single-trait matrices in memory and is therefore much faster and memory
efficient than gibbsf90 for multiple-tarit and correlated-effect models.
October 18
Number of missing records computed incorrectly in gibbs1f90.
January 11,2001
Assignment of hash to ija variable caused changes to the hash variable
declared as intent(in), resulting in error messages with some compilers.
During conversion, the hash variable is no longer modified.
Added option to addm(a,i,j,x,'f'), where x is the hash matrix, to store it in
full. Conversion of this matrix to the ija format is also in full.
January 22,2001
Added optional printing to a file in printmat in DENSEOP
January 24, 2001
Functions gen_wishart and gen_invwishart had nonstandard dispersion
argument: standarized and inverted. Now, the definittion of these functions
is compatible with literature. To obtain samples with quadratic form
ee and degrees of freedom df, use:
r=gen_invwishart(finverse_s(ee),df)
Applications using these functions: gibbsf90 and gibbs1f90, have been
modified accordingly.
Gibbs1f90 can store every n-th sample in file gibbs_samples.
Degrees of freedom for the iverted wishart generators in the gibbs
programs were slightly off in models with correlated effects
March 8, 2001
Gibbs1f90 did not handle correctly split fixed effects
Added gibbs2f90, which samples correlated effects jointly. This
results in much faster mixing in models with correlated random effects
such as maternal and random-regression.
March 12, 2001
IN REML and Gibbs programs, variance-covariances matrices with zeroed rows and
columns were found non semi-positive definite. This was dues to rounding errors
in a routine which calculates eigenvalues. Now, a matrix is declared
as non semi-positive definite when eigenvalues < -1e-15 rather than < 0.
March 22, 2001
REMLF90 was very slow when the number of effects was high, especially
with random regression models.
March 26, 2001
Small variances in REMLF90 were unnecessarily corrected for negative eigenvalues;
decreased 1d-4 in pos_def to 1d-10.
POS_DEF corrected to ignore zero rows/columns combinations.
May 11, 2001
POS-DEF did not compile with some compilers. Added keyword "recursive" in front
of it.
May 29, 2001
Gibbs family of programs did not work with Fujitsu compiler on Sun. This
compiler does not initialize dynamic structures by default. Added
explicit initializations to elements of ainv().
June 10, 2001
With some compilers (incl. Lahey and DEC), pedigree files were not read or
were overwritten if RANDOM_GROUP started with efefct 5 or 6. This was caused by
compiler bug/feature where unit 5 and 6 could not be open for file access.
File allocation has been changed so that units 5 and 6 are no longer used.
See blups1, the bottom of model.f90, and blupf90.f90 for details.
August 20,2001 (Changes made by Shogo Tsuruta)
In DENSEOP, cholesky decomposition, inversion and related subroutines/functions
were made faster by replacing dot_product with do loops. Speedup on matrices of 18x18 is over 2 times.
Gibbs2f90 speeded up by calculating the quadratic forms with IJA instead of
HASHM matrices
August 23, 2001
PRINTMAT printed text to standard output even if un was present.
August 30,2000
BLUPF90 modified to correct non-positive definite matrices of
variance components; negative eigenvalues are replaced by 1d-10.
November 7, 2001
Corrected formats for variances so that negative covariances are
not written as *******.
Decemeber 3, 2001
Modified POS_DEF subroutine so that the minimum eigenvalue is calculated relative to
the maximum eigenvalue, and all eigenvalues below a limit are corrected; see denseop
manual for reference. As a result, many programs and particularly AIREMLF90 and GIBBS*f90
are more stable when (co)variance matrices are close to singularity.
Makefile for Absoft verions 7.0 and later changed for higher optimization level and f95
comatibility.
January 10, 2002
BLUPF90 with the default PCG solver produced inaccurate solutions when variances
were large.
February 22, 2002
Parameter files for BLUPF90 and BLUP90IOD can now accept directives to change
default convergence criterion and teh maximum number of rounds; additionally,
parameter file for BLUPF90 can accept directives for solving methodology
and relaxation factor for SOR.
March 11
Gibbsf90 now writes samples to a file so that its output can be analyzed by
postgibbsf90.
May 16
cblup90reml now works more reliably for a single-trait threshold model.
function diag() now works also with single precision variables
June 14
AIREMLF90 now works with maternal and random regression models (fixes by
Tom Druet, INRA).
June 21
cblup90reml works faster with single-trait threshold models.
Sept 9
AIREMLF90 no longer generates
function diag() now works also with single precision variables
June 14
AIREMLF90 now works with maternal and random regression models (fixes by
Tom Druet, INRA).
June 21
cblup90reml works faster with single-trait threshold models.
Sept 9
AIREMLF90 no longer generates a message "wrong description of nested covariable'
Sept 20
Added Gibbs3f90, a version of gibbs2f90 with support for heterogeneous
residual variances.
Sept 24
In eigenvalue/eigenvector computations in denseop.f90, in rare instances
nonzero columns were recognized as 0 (Thanks Jan-Thijs van Kaam).
Oct 20
Nested covariables with nesting value of 0 do not work as missing value
(Detected by Tom Druet). Fix not yet applied. A workaround would be
to set covariables to 0 when nesting variable is also 0.
Nov 14
Renumf90 printed some statistics as *** or NaN. ALso some correlatiosn reported
as NaN.
Remlf90 printed correlations for nonestimated covariances as NaN.
Dec 9
In remlf90, some formats were too small; parameters for nesting
were incorrect in some multiple-trait cases.
Dec 12
RENUMF90 computed statistics using real numbers truncated to integers
Dec 13
Nested covariables with nesting value of 0 should work correctly
Feb 17, 2003
Pointer variables in sparse.f90 and fspak90.f90 nullified in declaration
Added support of user-defined relationships to remlf90
March 24
Missing values in gibbs3f90 were not correctly predicted.
(Thanks to Andres Legarra)
April 16
kron function moved from gibbs2 and gibbs3 to libs/denseop.f90
June 16
S. Tsuruta aded AIC to remlf90 and airemlf90
July 28
S. Tsuruta added options for convergence criteria to remlf90
July 31
Computations with dominance or user defined matrix may be inaccurate if these
effects are not last. Corrected in remlf90 only. Changes
include replacing neff with eff in add_g_dom, add_g_usr.
Oct 15
Read buffer in renumf90 was increased from 200 to 300 characters. Data lines
with greater length lead to long lines incorrectly read.