diff --git a/aprod.c b/aprod.c index 1e1507b99a8fa4a9dde4c0bf9e86345b453d4c1b..1b5051602a3d1dce65e88a4042c447b1026529ad 100644 Binary files a/aprod.c and b/aprod.c differ diff --git a/cblas.h b/cblas.h index c8233b444751846dfa1bc784824685ebfbef5ef0..2c282763fc7edd786295d20b8b16f5463243eac0 100644 --- a/cblas.h +++ b/cblas.h @@ -50,7 +50,7 @@ cblas_daxpy(const int N, const double alpha, const double *X, void cblas_dcopy(const long int N, const double *X, const int incX, - double *Y, const int incY); + double *Y, const int incY, int ntasks=0); double @@ -58,7 +58,7 @@ cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY); double -cblas_dnrm2(const long int N, const double *X, const int incX); +cblas_dnrm2(const long int N, const double *X, const int incX, int ntasks=0); void cblas_dscal(const long int N, const double alpha, double *X, const int incX); diff --git a/lsqr.c b/lsqr.c index 346591de1a0966c0aad93301ed5c0d690e02aa42..d1c6de0d228a52fbe1f6616ef0d66300ab6f67b0 100644 --- a/lsqr.c +++ b/lsqr.c @@ -32,10 +32,10 @@ Parallelized for ESA Gaia Mission. U. becciani A. Vecchiato 2013 // #include "cblas.h" //#endif -#define ZERO 0.0 -#define ONE 1.0 +#define ZERO 0.0 +#define ONE 1.0 void aprod(int mode, long int m, long int n, double x[], double y[], - double *ra,long int *matrixIndex, int *instrCol,int *instrConstrIlung, struct comData comlsqr,time_t *ompSec); + double *ra, long int *matrixIndex, int *instrCol, int *instrConstrIlung, struct comData comlsqr, time_t *ompSec); // --------------------------------------------------------------------- // d2norm returns sqrt( a**2 + b**2 ) with precautions // to avoid overflow. @@ -43,1328 +43,1476 @@ void aprod(int mode, long int m, long int n, double x[], double y[], // 21 Mar 1990: First version. // --------------------------------------------------------------------- static double -d2norm( const double a, const double b ) +d2norm(const double a, const double b) { double scale; const double zero = 0.0; - scale = fabs( a ) + fabs( b ); + scale = fabs(a) + fabs(b); if (scale == zero) return zero; else - return scale * sqrt( (a/scale)*(a/scale) + (b/scale)*(b/scale) ); + return scale * sqrt((a / scale) * (a / scale) + (b / scale) * (b / scale)); } static void -dload( const long int n, const double alpha, double x[] ) -{ +dload(const long int n, const double alpha, double x[]) +{ long int i; -#pragma omp for - for (i = 0; i < n; i++) x[i] = alpha; +/// FV_ EDIT + ///#pragma omp for + for (i = 0; i < n; i++) + x[i] = alpha; return; } // --------------------------------------------------------------------- // LSQR // --------------------------------------------------------------------- -void lsqr( - long int m, - long int n, -// void (*aprod)(int mode, int m, int n, double x[], double y[], -// void *UsrWrk), - double damp, -// void *UsrWrk, - double *knownTerms, // len = m reported as u - double *vVect, // len = n reported as v - double *wVect, // len = n reported as w - double *xSolution, // len = n reported as x - double *standardError, // len at least n. May be NULL. reported as se - double atol, - double btol, - double conlim, - int itnlim, - // The remaining variables are output only. - int *istop_out, - int *itn_out, - double *anorm_out, - double *acond_out, - double *rnorm_out, - double *arnorm_out, - double *xnorm_out, - double *systemMatrix, // reported as a - long int *matrixIndex, // reported as janew - int *instrCol, - int *instrConstrIlung, - double *preCondVect, - struct comData comlsqr) -{ -// ------------------------------------------------------------------ -// -// LSQR finds a solution x to the following problems: -// -// 1. Unsymmetric equations -- solve A*x = b -// -// 2. Linear least squares -- solve A*x = b -// in the least-squares sense -// -// 3. Damped least squares -- solve ( A )*x = ( b ) -// ( damp*I ) ( 0 ) -// in the least-squares sense -// -// where A is a matrix with m rows and n columns, b is an -// m-vector, and damp is a scalar. (All quantities are real.) -// The matrix A is intended to be large and sparse. It is accessed -// by means of subroutine calls of the form -// -// aprod ( mode, m, n, x, y, UsrWrk ) -// -// which must perform the following functions: -// -// If mode = 1, compute y = y + A*x. -// If mode = 2, compute x = x + A(transpose)*y. -// -// The vectors x and y are input parameters in both cases. -// If mode = 1, y should be altered without changing x. -// If mode = 2, x should be altered without changing y. -// The parameter UsrWrk may be used for workspace as described -// below. -// -// The rhs vector b is input via u, and subsequently overwritten. -// -// -// Note: LSQR uses an iterative method to approximate the solution. -// The number of iterations required to reach a certain accuracy -// depends strongly on the scaling of the problem. Poor scaling of -// the rows or columns of A should therefore be avoided where -// possible. -// -// For example, in problem 1 the solution is unaltered by -// row-scaling. If a row of A is very small or large compared to -// the other rows of A, the corresponding row of ( A b ) should be -// scaled up or down. -// -// In problems 1 and 2, the solution x is easily recovered -// following column-scaling. Unless better information is known, -// the nonzero columns of A should be scaled so that they all have -// the same Euclidean norm (e.g., 1.0). -// -// In problem 3, there is no freedom to re-scale if damp is -// nonzero. However, the value of damp should be assigned only -// after attention has been paid to the scaling of A. -// -// The parameter damp is intended to help regularize -// ill-conditioned systems, by preventing the true solution from -// being very large. Another aid to regularization is provided by -// the parameter acond, which may be used to terminate iterations -// before the computed solution becomes very large. -// -// Note that x is not an input parameter. -// If some initial estimate x0 is known and if damp = 0, -// one could proceed as follows: -// -// 1. Compute a residual vector r0 = b - A*x0. -// 2. Use LSQR to solve the system A*dx = r0. -// 3. Add the correction dx to obtain a final solution x = x0 + dx. -// -// This requires that x0 be available before and after the call -// to LSQR. To judge the benefits, suppose LSQR takes k1 iterations -// to solve A*x = b and k2 iterations to solve A*dx = r0. -// If x0 is "good", norm(r0) will be smaller than norm(b). -// If the same stopping tolerances atol and btol are used for each -// system, k1 and k2 will be similar, but the final solution x0 + dx -// should be more accurate. The only way to reduce the total work -// is to use a larger stopping tolerance for the second system. -// If some value btol is suitable for A*x = b, the larger value -// btol*norm(b)/norm(r0) should be suitable for A*dx = r0. -// -// Preconditioning is another way to reduce the number of iterations. -// If it is possible to solve a related system M*x = b efficiently, -// where M approximates A in some helpful way -// (e.g. M - A has low rank or its elements are small relative to -// those of A), LSQR may converge more rapidly on the system -// A*M(inverse)*z = b, -// after which x can be recovered by solving M*x = z. -// -// NOTE: If A is symmetric, LSQR should not be used! -// Alternatives are the symmetric conjugate-gradient method (cg) -// and/or SYMMLQ. -// SYMMLQ is an implementation of symmetric cg that applies to -// any symmetric A and will converge more rapidly than LSQR. -// If A is positive definite, there are other implementations of -// symmetric cg that require slightly less work per iteration -// than SYMMLQ (but will take the same number of iterations). -// -// -// Notation -// -------- -// -// The following quantities are used in discussing the subroutine -// parameters: -// -// Abar = ( A ), bbar = ( b ) -// ( damp*I ) ( 0 ) -// -// r = b - A*x, rbar = bbar - Abar*x -// -// rnorm = sqrt( norm(r)**2 + damp**2 * norm(x)**2 ) -// = norm( rbar ) -// -// relpr = the relative precision of floating-point arithmetic -// on the machine being used. On most machines, -// relpr is about 1.0e-7 and 1.0d-16 in single and double -// precision respectively. -// -// LSQR minimizes the function rnorm with respect to x. -// -// -// Parameters -// ---------- -// -// m input m, the number of rows in A. -// -// n input n, the number of columns in A. -// -// aprod external See above. -// -// damp input The damping parameter for problem 3 above. -// (damp should be 0.0 for problems 1 and 2.) -// If the system A*x = b is incompatible, values -// of damp in the range 0 to sqrt(relpr)*norm(A) -// will probably have a negligible effect. -// Larger values of damp will tend to decrease -// the norm of x and reduce the number of -// iterations required by LSQR. -// -// The work per iteration and the storage needed -// by LSQR are the same for all values of damp. -// -// rw workspace Transit pointer to user's workspace. -// Note: LSQR does not explicitly use this -// parameter, but passes it to subroutine aprod for -// possible use as workspace. -// -// u(m) input The rhs vector b. Beware that u is -// over-written by LSQR. -// -// v(n) workspace -// -// w(n) workspace -// -// x(n) output Returns the computed solution x. -// -// se(*) output If m .gt. n or damp .gt. 0, the system is -// (maybe) overdetermined and the standard errors may be -// useful. (See the first LSQR reference.) -// Otherwise (m .le. n and damp = 0) they do not -// mean much. Some time and storage can be saved -// by setting se = NULL. In that case, se will -// not be touched. -// -// If se is not NULL, then the dimension of se must -// be n or more. se(1:n) then returns standard error -// estimates for the components of x. -// For each i, se(i) is set to the value -// rnorm * sqrt( sigma(i,i) / t ), -// where sigma(i,i) is an estimate of the i-th -// diagonal of the inverse of Abar(transpose)*Abar -// and t = 1 if m .le. n, -// t = m - n if m .gt. n and damp = 0, -// t = m if damp .ne. 0. -// -// atol input An estimate of the relative error in the data -// defining the matrix A. For example, -// if A is accurate to about 6 digits, set -// atol = 1.0e-6 . -// -// btol input An estimate of the relative error in the data -// defining the rhs vector b. For example, -// if b is accurate to about 6 digits, set -// btol = 1.0e-6 . -// -// conlim input An upper limit on cond(Abar), the apparent -// condition number of the matrix Abar. -// Iterations will be terminated if a computed -// estimate of cond(Abar) exceeds conlim. -// This is intended to prevent certain small or -// zero singular values of A or Abar from -// coming into effect and causing unwanted growth -// in the computed solution. -// -// conlim and damp may be used separately or -// together to regularize ill-conditioned systems. -// -// Normally, conlim should be in the range -// 1000 to 1/relpr. -// Suggested value: -// conlim = 1/(100*relpr) for compatible systems, -// conlim = 1/(10*sqrt(relpr)) for least squares. -// -// Note: If the user is not concerned about the parameters -// atol, btol and conlim, any or all of them may be set -// to zero. The effect will be the same as the values -// relpr, relpr and 1/relpr respectively. -// -// itnlim input An upper limit on the number of iterations. -// Suggested value: -// itnlim = n/2 for well-conditioned systems -// with clustered singular values, -// itnlim = 4*n otherwise. -// -// nout input File number for printed output. If positive, -// a summary will be printed on file nout. -// -// istop output An integer giving the reason for termination: -// -// 0 x = 0 is the exact solution. -// No iterations were performed. -// -// 1 The equations A*x = b are probably -// compatible. Norm(A*x - b) is sufficiently -// small, given the values of atol and btol. -// -// 2 damp is zero. The system A*x = b is probably -// not compatible. A least-squares solution has -// been obtained that is sufficiently accurate, -// given the value of atol. -// -// 3 damp is nonzero. A damped least-squares -// solution has been obtained that is sufficiently -// accurate, given the value of atol. -// -// 4 An estimate of cond(Abar) has exceeded -// conlim. The system A*x = b appears to be -// ill-conditioned. Otherwise, there could be an -// error in subroutine aprod. -// -// 5 The iteration limit itnlim was reached. -// -// itn output The number of iterations performed. -// -// anorm output An estimate of the Frobenius norm of Abar. -// This is the square-root of the sum of squares -// of the elements of Abar. -// If damp is small and if the columns of A -// have all been scaled to have length 1.0, -// anorm should increase to roughly sqrt(n). -// A radically different value for anorm may -// indicate an error in subroutine aprod (there -// may be an inconsistency between modes 1 and 2). -// -// acond output An estimate of cond(Abar), the condition -// number of Abar. A very high value of acond -// may again indicate an error in aprod. -// -// rnorm output An estimate of the final value of norm(rbar), -// the function being minimized (see notation -// above). This will be small if A*x = b has -// a solution. -// -// arnorm output An estimate of the final value of -// norm( Abar(transpose)*rbar ), the norm of -// the residual for the usual normal equations. -// This should be small in all cases. (arnorm -// will often be smaller than the true value -// computed from the output vector x.) -// -// xnorm output An estimate of the norm of the final -// solution vector x. -// -// -// Subroutines and functions used -// ------------------------------ -// -// USER aprod -// CBLAS dcopy, dnrm2, dscal (see Lawson et al. below) -// -// -// References -// ---------- -// -// C.C. Paige and M.A. Saunders, LSQR: An algorithm for sparse -// linear equations and sparse least squares, -// ACM Transactions on Mathematical Software 8, 1 (March 1982), -// pp. 43-71. -// -// C.C. Paige and M.A. Saunders, Algorithm 583, LSQR: Sparse -// linear equations and least-squares problems, -// ACM Transactions on Mathematical Software 8, 2 (June 1982), -// pp. 195-209. -// -// C.L. Lawson, R.J. Hanson, D.R. Kincaid and F.T. Krogh, -// Basic linear algebra subprograms for Fortran usage, -// ACM Transactions on Mathematical Software 5, 3 (Sept 1979), -// pp. 308-323 and 324-325. -// ------------------------------------------------------------------ -// -// -// LSQR development: -// 22 Feb 1982: LSQR sent to ACM TOMS to become Algorithm 583. -// 15 Sep 1985: Final F66 version. LSQR sent to "misc" in netlib. -// 13 Oct 1987: Bug (Robert Davies, DSIR). Have to delete -// if ( (one + dabs(t)) .le. one ) GO TO 200 -// from loop 200. The test was an attempt to reduce -// underflows, but caused w(i) not to be updated. -// 17 Mar 1989: First F77 version. -// 04 May 1989: Bug (David Gay, AT&T). When the second beta is zero, -// rnorm = 0 and -// test2 = arnorm / (anorm * rnorm) overflows. -// Fixed by testing for rnorm = 0. -// 05 May 1989: Sent to "misc" in netlib. -// 14 Mar 1990: Bug (John Tomlin via IBM OSL testing). -// Setting rhbar2 = rhobar**2 + dampsq can give zero -// if rhobar underflows and damp = 0. -// Fixed by testing for damp = 0 specially. -// 15 Mar 1990: Converted to lower case. -// 21 Mar 1990: d2norm introduced to avoid overflow in numerous -// items like c = sqrt( a**2 + b**2 ). -// 04 Sep 1991: wantse added as an argument to LSQR, to make -// standard errors optional. This saves storage and -// time when se(*) is not wanted. -// 13 Feb 1992: istop now returns a value in [1,5], not [1,7]. -// 1, 2 or 3 means that x solves one of the problems -// Ax = b, min norm(Ax - b) or damped least squares. -// 4 means the limit on cond(A) was reached. -// 5 means the limit on iterations was reached. -// 07 Dec 1994: Keep track of dxmax = max_k norm( phi_k * d_k ). -// So far, this is just printed at the end. -// A large value (relative to norm(x)) indicates -// significant cancellation in forming -// x = D*f = sum( phi_k * d_k ). -// A large column of D need NOT be serious if the -// corresponding phi_k is small. -// 27 Dec 1994: Include estimate of alfa_opt in iteration log. -// alfa_opt is the optimal scale factor for the -// residual in the "augmented system", as described by -// A. Bjorck (1992), -// Pivoting and stability in the augmented system method, -// in D. F. Griffiths and G. A. Watson (eds.), -// "Numerical Analysis 1991", -// Proceedings of the 14th Dundee Conference, -// Pitman Research Notes in Mathematics 260, -// Longman Scientific and Technical, Harlow, Essex, 1992. -// 14 Apr 2006: "Line-by-line" conversion to ISO C by -// Michael P. Friedlander. -// -// -// Michael A. Saunders mike@sol-michael.stanford.edu -// Dept of Operations Research na.Msaunders@na-net.ornl.gov -// Stanford University -// Stanford, CA 94305-4022 (415) 723-1875 -//----------------------------------------------------------------------- +void lsqr( + long int m, + long int n, + // void (*aprod)(int mode, int m, int n, double x[], double y[], + // void *UsrWrk), + double damp, + // void *UsrWrk, + double *knownTerms, // len = m reported as u + double *vVect, // len = n reported as v + double *wVect, // len = n reported as w + double *xSolution, // len = n reported as x + double *standardError, // len at least n. May be NULL. reported as se + double atol, + double btol, + double conlim, + int itnlim, + // The remaining variables are output only. + int *istop_out, + int *itn_out, + double *anorm_out, + double *acond_out, + double *rnorm_out, + double *arnorm_out, + double *xnorm_out, + double *systemMatrix, // reported as a + long int *matrixIndex, // reported as janew + int *instrCol, + int *instrConstrIlung, + double *preCondVect, + struct comData comlsqr) + { + // ------------------------------------------------------------------ + // + // LSQR finds a solution x to the following problems: + // + // 1. Unsymmetric equations -- solve A*x = b + // + // 2. Linear least squares -- solve A*x = b + // in the least-squares sense + // + // 3. Damped least squares -- solve ( A )*x = ( b ) + // ( damp*I ) ( 0 ) + // in the least-squares sense + // + // where A is a matrix with m rows and n columns, b is an + // m-vector, and damp is a scalar. (All quantities are real.) + // The matrix A is intended to be large and sparse. It is accessed + // by means of subroutine calls of the form + // + // aprod ( mode, m, n, x, y, UsrWrk ) + // + // which must perform the following functions: + // + // If mode = 1, compute y = y + A*x. + // If mode = 2, compute x = x + A(transpose)*y. + // + // The vectors x and y are input parameters in both cases. + // If mode = 1, y should be altered without changing x. + // If mode = 2, x should be altered without changing y. + // The parameter UsrWrk may be used for workspace as described + // below. + // + // The rhs vector b is input via u, and subsequently overwritten. + // + // + // Note: LSQR uses an iterative method to approximate the solution. + // The number of iterations required to reach a certain accuracy + // depends strongly on the scaling of the problem. Poor scaling of + // the rows or columns of A should therefore be avoided where + // possible. + // + // For example, in problem 1 the solution is unaltered by + // row-scaling. If a row of A is very small or large compared to + // the other rows of A, the corresponding row of ( A b ) should be + // scaled up or down. + // + // In problems 1 and 2, the solution x is easily recovered + // following column-scaling. Unless better information is known, + // the nonzero columns of A should be scaled so that they all have + // the same Euclidean norm (e.g., 1.0). + // + // In problem 3, there is no freedom to re-scale if damp is + // nonzero. However, the value of damp should be assigned only + // after attention has been paid to the scaling of A. + // + // The parameter damp is intended to help regularize + // ill-conditioned systems, by preventing the true solution from + // being very large. Another aid to regularization is provided by + // the parameter acond, which may be used to terminate iterations + // before the computed solution becomes very large. + // + // Note that x is not an input parameter. + // If some initial estimate x0 is known and if damp = 0, + // one could proceed as follows: + // + // 1. Compute a residual vector r0 = b - A*x0. + // 2. Use LSQR to solve the system A*dx = r0. + // 3. Add the correction dx to obtain a final solution x = x0 + dx. + // + // This requires that x0 be available before and after the call + // to LSQR. To judge the benefits, suppose LSQR takes k1 iterations + // to solve A*x = b and k2 iterations to solve A*dx = r0. + // If x0 is "good", norm(r0) will be smaller than norm(b). + // If the same stopping tolerances atol and btol are used for each + // system, k1 and k2 will be similar, but the final solution x0 + dx + // should be more accurate. The only way to reduce the total work + // is to use a larger stopping tolerance for the second system. + // If some value btol is suitable for A*x = b, the larger value + // btol*norm(b)/norm(r0) should be suitable for A*dx = r0. + // + // Preconditioning is another way to reduce the number of iterations. + // If it is possible to solve a related system M*x = b efficiently, + // where M approximates A in some helpful way + // (e.g. M - A has low rank or its elements are small relative to + // those of A), LSQR may converge more rapidly on the system + // A*M(inverse)*z = b, + // after which x can be recovered by solving M*x = z. + // + // NOTE: If A is symmetric, LSQR should not be used! + // Alternatives are the symmetric conjugate-gradient method (cg) + // and/or SYMMLQ. + // SYMMLQ is an implementation of symmetric cg that applies to + // any symmetric A and will converge more rapidly than LSQR. + // If A is positive definite, there are other implementations of + // symmetric cg that require slightly less work per iteration + // than SYMMLQ (but will take the same number of iterations). + // + // + // Notation + // -------- + // + // The following quantities are used in discussing the subroutine + // parameters: + // + // Abar = ( A ), bbar = ( b ) + // ( damp*I ) ( 0 ) + // + // r = b - A*x, rbar = bbar - Abar*x + // + // rnorm = sqrt( norm(r)**2 + damp**2 * norm(x)**2 ) + // = norm( rbar ) + // + // relpr = the relative precision of floating-point arithmetic + // on the machine being used. On most machines, + // relpr is about 1.0e-7 and 1.0d-16 in single and double + // precision respectively. + // + // LSQR minimizes the function rnorm with respect to x. + // + // + // Parameters + // ---------- + // + // m input m, the number of rows in A. + // + // n input n, the number of columns in A. + // + // aprod external See above. + // + // damp input The damping parameter for problem 3 above. + // (damp should be 0.0 for problems 1 and 2.) + // If the system A*x = b is incompatible, values + // of damp in the range 0 to sqrt(relpr)*norm(A) + // will probably have a negligible effect. + // Larger values of damp will tend to decrease + // the norm of x and reduce the number of + // iterations required by LSQR. + // + // The work per iteration and the storage needed + // by LSQR are the same for all values of damp. + // + // rw workspace Transit pointer to user's workspace. + // Note: LSQR does not explicitly use this + // parameter, but passes it to subroutine aprod for + // possible use as workspace. + // + // u(m) input The rhs vector b. Beware that u is + // over-written by LSQR. + // + // v(n) workspace + // + // w(n) workspace + // + // x(n) output Returns the computed solution x. + // + // se(*) output If m .gt. n or damp .gt. 0, the system is + // (maybe) overdetermined and the standard errors may be + // useful. (See the first LSQR reference.) + // Otherwise (m .le. n and damp = 0) they do not + // mean much. Some time and storage can be saved + // by setting se = NULL. In that case, se will + // not be touched. + // + // If se is not NULL, then the dimension of se must + // be n or more. se(1:n) then returns standard error + // estimates for the components of x. + // For each i, se(i) is set to the value + // rnorm * sqrt( sigma(i,i) / t ), + // where sigma(i,i) is an estimate of the i-th + // diagonal of the inverse of Abar(transpose)*Abar + // and t = 1 if m .le. n, + // t = m - n if m .gt. n and damp = 0, + // t = m if damp .ne. 0. + // + // atol input An estimate of the relative error in the data + // defining the matrix A. For example, + // if A is accurate to about 6 digits, set + // atol = 1.0e-6 . + // + // btol input An estimate of the relative error in the data + // defining the rhs vector b. For example, + // if b is accurate to about 6 digits, set + // btol = 1.0e-6 . + // + // conlim input An upper limit on cond(Abar), the apparent + // condition number of the matrix Abar. + // Iterations will be terminated if a computed + // estimate of cond(Abar) exceeds conlim. + // This is intended to prevent certain small or + // zero singular values of A or Abar from + // coming into effect and causing unwanted growth + // in the computed solution. + // + // conlim and damp may be used separately or + // together to regularize ill-conditioned systems. + // + // Normally, conlim should be in the range + // 1000 to 1/relpr. + // Suggested value: + // conlim = 1/(100*relpr) for compatible systems, + // conlim = 1/(10*sqrt(relpr)) for least squares. + // + // Note: If the user is not concerned about the parameters + // atol, btol and conlim, any or all of them may be set + // to zero. The effect will be the same as the values + // relpr, relpr and 1/relpr respectively. + // + // itnlim input An upper limit on the number of iterations. + // Suggested value: + // itnlim = n/2 for well-conditioned systems + // with clustered singular values, + // itnlim = 4*n otherwise. + // + // nout input File number for printed output. If positive, + // a summary will be printed on file nout. + // + // istop output An integer giving the reason for termination: + // + // 0 x = 0 is the exact solution. + // No iterations were performed. + // + // 1 The equations A*x = b are probably + // compatible. Norm(A*x - b) is sufficiently + // small, given the values of atol and btol. + // + // 2 damp is zero. The system A*x = b is probably + // not compatible. A least-squares solution has + // been obtained that is sufficiently accurate, + // given the value of atol. + // + // 3 damp is nonzero. A damped least-squares + // solution has been obtained that is sufficiently + // accurate, given the value of atol. + // + // 4 An estimate of cond(Abar) has exceeded + // conlim. The system A*x = b appears to be + // ill-conditioned. Otherwise, there could be an + // error in subroutine aprod. + // + // 5 The iteration limit itnlim was reached. + // + // itn output The number of iterations performed. + // + // anorm output An estimate of the Frobenius norm of Abar. + // This is the square-root of the sum of squares + // of the elements of Abar. + // If damp is small and if the columns of A + // have all been scaled to have length 1.0, + // anorm should increase to roughly sqrt(n). + // A radically different value for anorm may + // indicate an error in subroutine aprod (there + // may be an inconsistency between modes 1 and 2). + // + // acond output An estimate of cond(Abar), the condition + // number of Abar. A very high value of acond + // may again indicate an error in aprod. + // + // rnorm output An estimate of the final value of norm(rbar), + // the function being minimized (see notation + // above). This will be small if A*x = b has + // a solution. + // + // arnorm output An estimate of the final value of + // norm( Abar(transpose)*rbar ), the norm of + // the residual for the usual normal equations. + // This should be small in all cases. (arnorm + // will often be smaller than the true value + // computed from the output vector x.) + // + // xnorm output An estimate of the norm of the final + // solution vector x. + // + // + // Subroutines and functions used + // ------------------------------ + // + // USER aprod + // CBLAS dcopy, dnrm2, dscal (see Lawson et al. below) + // + // + // References + // ---------- + // + // C.C. Paige and M.A. Saunders, LSQR: An algorithm for sparse + // linear equations and sparse least squares, + // ACM Transactions on Mathematical Software 8, 1 (March 1982), + // pp. 43-71. + // + // C.C. Paige and M.A. Saunders, Algorithm 583, LSQR: Sparse + // linear equations and least-squares problems, + // ACM Transactions on Mathematical Software 8, 2 (June 1982), + // pp. 195-209. + // + // C.L. Lawson, R.J. Hanson, D.R. Kincaid and F.T. Krogh, + // Basic linear algebra subprograms for Fortran usage, + // ACM Transactions on Mathematical Software 5, 3 (Sept 1979), + // pp. 308-323 and 324-325. + // ------------------------------------------------------------------ + // + // + // LSQR development: + // 22 Feb 1982: LSQR sent to ACM TOMS to become Algorithm 583. + // 15 Sep 1985: Final F66 version. LSQR sent to "misc" in netlib. + // 13 Oct 1987: Bug (Robert Davies, DSIR). Have to delete + // if ( (one + dabs(t)) .le. one ) GO TO 200 + // from loop 200. The test was an attempt to reduce + // underflows, but caused w(i) not to be updated. + // 17 Mar 1989: First F77 version. + // 04 May 1989: Bug (David Gay, AT&T). When the second beta is zero, + // rnorm = 0 and + // test2 = arnorm / (anorm * rnorm) overflows. + // Fixed by testing for rnorm = 0. + // 05 May 1989: Sent to "misc" in netlib. + // 14 Mar 1990: Bug (John Tomlin via IBM OSL testing). + // Setting rhbar2 = rhobar**2 + dampsq can give zero + // if rhobar underflows and damp = 0. + // Fixed by testing for damp = 0 specially. + // 15 Mar 1990: Converted to lower case. + // 21 Mar 1990: d2norm introduced to avoid overflow in numerous + // items like c = sqrt( a**2 + b**2 ). + // 04 Sep 1991: wantse added as an argument to LSQR, to make + // standard errors optional. This saves storage and + // time when se(*) is not wanted. + // 13 Feb 1992: istop now returns a value in [1,5], not [1,7]. + // 1, 2 or 3 means that x solves one of the problems + // Ax = b, min norm(Ax - b) or damped least squares. + // 4 means the limit on cond(A) was reached. + // 5 means the limit on iterations was reached. + // 07 Dec 1994: Keep track of dxmax = max_k norm( phi_k * d_k ). + // So far, this is just printed at the end. + // A large value (relative to norm(x)) indicates + // significant cancellation in forming + // x = D*f = sum( phi_k * d_k ). + // A large column of D need NOT be serious if the + // corresponding phi_k is small. + // 27 Dec 1994: Include estimate of alfa_opt in iteration log. + // alfa_opt is the optimal scale factor for the + // residual in the "augmented system", as described by + // A. Bjorck (1992), + // Pivoting and stability in the augmented system method, + // in D. F. Griffiths and G. A. Watson (eds.), + // "Numerical Analysis 1991", + // Proceedings of the 14th Dundee Conference, + // Pitman Research Notes in Mathematics 260, + // Longman Scientific and Technical, Harlow, Essex, 1992. + // 14 Apr 2006: "Line-by-line" conversion to ISO C by + // Michael P. Friedlander. + // + // + // Michael A. Saunders mike@sol-michael.stanford.edu + // Dept of Operations Research na.Msaunders@na-net.ornl.gov + // Stanford University + // Stanford, CA 94305-4022 (415) 723-1875 + //----------------------------------------------------------------------- + + // Local copies of output variables. Output vars are assigned at exit. -// Local copies of output variables. Output vars are assigned at exit. + + + - int - istop = 0, - itn = 0; - double - anorm = ZERO, - acond = ZERO, - rnorm = ZERO, - arnorm = ZERO, - xnorm = ZERO; - + + + int istop = 0; + int itn = 0; + + +double + anorm = ZERO, + acond = ZERO, + rnorm = ZERO, + arnorm = ZERO, + xnorm = ZERO; + // Local variables - const bool - extra = false, // true for extra printing below. - damped = damp > ZERO, - wantse = standardError != NULL; - long int i; - int - maxdx, nconv, nstop; - double - alfopt, alpha, arnorm0, beta, bnorm, - cs, cs1, cs2, ctol, - delta, dknorm, dknormSum, dnorm, dxk, dxmax, - gamma, gambar, phi, phibar, psi, - res2, rho, rhobar, rhbar1, - rhs, rtol, sn, sn1, sn2, - t, tau, temp, test1, test2, test3, - theta, t1, t2, t3, xnorm1, z, zbar; - char - enter[] = "Enter LSQR. ", - exitLsqr[] = "Exit LSQR. ", - msg[6][100] = +const bool + extra = false, // true for extra printing below. + damped = damp > ZERO, + wantse = standardError != NULL; +long int i; +int + maxdx, + nconv, nstop; +double + alfopt, + alpha, arnorm0, beta, bnorm, + cs, cs1, cs2, ctol, + delta, dknorm, dknormSum, dnorm, dxk, dxmax, + gamma, gambar, phi, phibar, psi, + res2, rho, rhobar, rhbar1, + rhs, rtol, sn, sn1, sn2, + t, tau, temp, test1, test2, test3, + theta, t1, t2, t3, xnorm1, z, zbar; +char + enter[] = "Enter LSQR. ", + exitLsqr[] = "Exit LSQR. ", + msg[6][100] = { {"The exact solution is x = 0"}, {"A solution to Ax = b was found, given atol, btol"}, {"A least-squares solution was found, given atol"}, {"A damped least-squares solution was found, given atol"}, {"Cond(Abar) seems to be too large, given conlim"}, - {"The iteration limit was reached"} - }; - char lsqrOut[] = "lsqr-output"; // (buffer flush problem) - char wpath[1024]; - size_t sizePath=1020; - + {"The iteration limit was reached"}}; +char lsqrOut[] = "lsqr-output"; // (buffer flush problem) +char wpath[1024]; +size_t sizePath = 1020; + //----------------------------------------------------------------------- // Format strings. - char fmt_1000[] = - " %s Least-squares solution of Ax = b\n" - " The matrix A has %ld rows and %ld columns\n" - " damp = %-22.2e wantse = %10i\n" - " atol = %-22.2e conlim = %10.2e\n" - " btol = %-22.2e itnlim = %10d\n\n"; - char fmt_1200[] = - " Itn x(1) Function" - " Compatible LS Norm A Cond A\n"; - char fmt_1300[] = - " Itn x(1) Function" - " Compatible LS Norm Abar Cond Abar\n"; - char fmt_1400[] = - " phi dknorm dxk alfa_opt\n"; - char fmt_1500_extra[] = - " %6d %16.9e %16.9e %9.2e %9.2e %8.1e %8.1e %8.1e %7.1e %7.1e %7.1e\n"; - char fmt_1500[] = - " %6d %16.9e %16.9e %9.2e %9.2e %8.1e %8.1e\n"; - char fmt_1550[] = - " %6d %16.9e %16.9e %9.2e %9.2e\n"; - char fmt_1600[] = - "\n"; - char fmt_2000[] = - "\n" - " %s istop = %-10d itn = %-10d\n" - " %s anorm = %11.5e acond = %11.5e\n" - " %s vnorm = %11.5e xnorm = %11.5e\n" - " %s rnorm = %11.5e arnorm = %11.5e\n"; - char fmt_2100[] = - " %s max dx = %7.1e occured at itn %-9d\n" - " %s = %7.1e*xnorm\n"; - char fmt_3000[] = - " %s %s\n"; +char fmt_1000[] = + " %s Least-squares solution of Ax = b\n" + " The matrix A has %ld rows and %ld columns\n" + " damp = %-22.2e wantse = %10i\n" + " atol = %-22.2e conlim = %10.2e\n" + " btol = %-22.2e itnlim = %10d\n\n"; +char fmt_1200[] = + " Itn x(1) Function" + " Compatible LS Norm A Cond A\n"; +char fmt_1300[] = + " Itn x(1) Function" + " Compatible LS Norm Abar Cond Abar\n"; +char fmt_1400[] = + " phi dknorm dxk alfa_opt\n"; +char fmt_1500_extra[] = + " %6d %16.9e %16.9e %9.2e %9.2e %8.1e %8.1e %8.1e %7.1e %7.1e %7.1e\n"; +char fmt_1500[] = + " %6d %16.9e %16.9e %9.2e %9.2e %8.1e %8.1e\n"; +char fmt_1550[] = + " %6d %16.9e %16.9e %9.2e %9.2e\n"; +char fmt_1600[] = + "\n"; +char fmt_2000[] = + "\n" + " %s istop = %-10d itn = %-10d\n" + " %s anorm = %11.5e acond = %11.5e\n" + " %s vnorm = %11.5e xnorm = %11.5e\n" + " %s rnorm = %11.5e arnorm = %11.5e\n"; +char fmt_2100[] = + " %s max dx = %7.1e occured at itn %-9d\n" + " %s = %7.1e*xnorm\n"; +char fmt_3000[] = + " %s %s\n"; ///////////// Specific definitions - time_t startCycleTime,endCycleTime,totTime,partialTime,CPRtimeStart,CPRtimeEnd,ompSec=0; - int writeCPR; //=0 no write CPR, =1 write CPR and continue, =2 write CPR and return - int noCPR; - int myid,nproc; - long int *mapNoss, *mapNcoeff; - MPI_Status status; - long int nunkSplit, localAstro, localAstroMax, other; - int nAstroElements; - int debugMode; - long VrIdAstroPDim, VrIdAstroPDimMax; - long nDegFreedomAtt; - int nAttAxes; - int nAttParam,nInstrParam,nGlobalParam,itnTest,nAttP,numOfExtStar,numOfBarStar,numOfExtAttCol,nobs; - long mapNossBefore; - short nAstroPSolved,nInstrPsolved; - double cycleStartMpiTime, cycleEndMpiTime; - double *vAuxVect; - int CPRCount; - FILE *fpCPRend, *nout, *fpItnLimitNew; - int nEqExtConstr,nEqBarConstr; - int nElemIC,nOfInstrConstr; - double *kAuxcopy; - double *kcopy; -//////////////////////////////// +time_t startCycleTime, endCycleTime, totTime, partialTime, CPRtimeStart, CPRtimeEnd, ompSec = 0; +int writeCPR; //=0 no write CPR, =1 write CPR and continue, =2 write CPR and return +int noCPR; +int myid, nproc; +long int *mapNoss, *mapNcoeff; +MPI_Status status; +long int nunkSplit, localAstro, localAstroMax, other; +int nAstroElements; +int debugMode; +long VrIdAstroPDim, VrIdAstroPDimMax; +long nDegFreedomAtt; +int nAttAxes; +int nAttParam, nInstrParam, nGlobalParam, itnTest, nAttP, numOfExtStar, numOfBarStar, numOfExtAttCol, nobs; +long mapNossBefore; +short nAstroPSolved, nInstrPsolved; +double cycleStartMpiTime, cycleEndMpiTime; +double *vAuxVect; +int CPRCount; +FILE *fpCPRend, *nout, *fpItnLimitNew; +int nEqExtConstr, nEqBarConstr; +int nElemIC, nOfInstrConstr; +double *kAuxcopy; +double *kcopy; +int ntasks; +int nthreads; +//////////////////////////////// // Initialize. -myid=comlsqr.myid; -nproc=comlsqr.nproc; -mapNcoeff=comlsqr.mapNcoeff; -mapNoss=comlsqr.mapNoss; -nAttParam=comlsqr.nAttParam; -nInstrParam=comlsqr.nInstrParam; -nGlobalParam=comlsqr.nGlobalParam; -nunkSplit=comlsqr.nunkSplit; -VrIdAstroPDim=comlsqr.VrIdAstroPDim; -VrIdAstroPDimMax=comlsqr.VrIdAstroPDimMax; -nAstroPSolved=comlsqr.nAstroPSolved; -nEqExtConstr=comlsqr.nEqExtConstr; -nEqBarConstr=comlsqr.nEqBarConstr; -nDegFreedomAtt=comlsqr.nDegFreedomAtt; -nAttAxes=comlsqr.nAttAxes; -localAstro= VrIdAstroPDim*nAstroPSolved; -localAstroMax=VrIdAstroPDimMax*nAstroPSolved; -numOfExtStar=comlsqr.numOfExtStar; -numOfBarStar=comlsqr.numOfBarStar; -nAttP=comlsqr.nAttP; -numOfExtAttCol=comlsqr.numOfExtAttCol; -mapNossBefore=comlsqr.mapNossBefore; -nobs=comlsqr.nobs; -debugMode=comlsqr.debugMode; -noCPR=comlsqr.noCPR; -nElemIC=comlsqr.nElemIC; -nOfInstrConstr=comlsqr.nOfInstrConstr; -nInstrPsolved=comlsqr.nInstrPSolved; - -other=(long)nAttParam + nInstrParam + nGlobalParam; -if(nAstroPSolved) vAuxVect=(double *) calloc(localAstroMax,sizeof(double)); - -comlsqr.itn=itn; -totTime=comlsqr.totSec; -partialTime=comlsqr.totSec; - - if(debugMode){ - for(long i=0; i ZERO) ctol = ONE / conlim; -anorm = ZERO; -acond = ZERO; -dnorm = ZERO; -dxmax = ZERO; -res2 = ZERO; -psi = ZERO; -xnorm = ZERO; -xnorm1 = ZERO; -cs2 = - ONE; -sn2 = ZERO; -z = ZERO; - -CPRCount=0; +itn = 0; +istop = 0; +nstop = 0; +maxdx = 0; +ctol = ZERO; +if (conlim > ZERO) + ctol = ONE / conlim; +anorm = ZERO; +acond = ZERO; +dnorm = ZERO; +dxmax = ZERO; +res2 = ZERO; +psi = ZERO; +xnorm = ZERO; +xnorm1 = ZERO; +cs2 = -ONE; +sn2 = ZERO; +z = ZERO; + + + + +CPRCount = 0; // ------------------------------------------------------------------ // Set up the first vectors u and v for the bidiagonalization. // These satisfy beta*u = b, alpha*v = A(transpose)*u. // ------------------------------------------------------------------ -dload( nunkSplit, 0.0, vVect ); -dload( nunkSplit, 0.0, xSolution ); -if ( wantse ) - dload( nunkSplit, 0.0, standardError ); -alpha = ZERO; +dload(nunkSplit, 0.0, vVect); +dload(nunkSplit, 0.0, xSolution); +if (wantse) + dload(nunkSplit, 0.0, standardError); +alpha = ZERO; // Find the maximum value on u - double betaLoc; - if(myid==0) betaLoc=cblas_dnrm2(mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr,knownTerms,1); - else betaLoc=cblas_dnrm2(mapNoss[myid],knownTerms,1); +double betaLoc; +if (myid == 0) + betaLoc = cblas_dnrm2(mapNoss[myid] + nEqExtConstr + nEqBarConstr + nOfInstrConstr, knownTerms, 1); +else + betaLoc = cblas_dnrm2(mapNoss[myid], knownTerms, 1); -double betaLoc2=betaLoc*betaLoc; +double betaLoc2 = betaLoc * betaLoc; double betaGlob; -MPI_Allreduce(&betaLoc2, &betaGlob,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); -beta=sqrt(betaGlob); -if (beta > ZERO) +MPI_Allreduce(&betaLoc2, &betaGlob, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +beta = sqrt(betaGlob); +if (beta > ZERO) { - cblas_dscal ( mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr, (ONE / beta), knownTerms, 1 ); - MPI_Barrier(MPI_COMM_WORLD); + cblas_dscal(mapNoss[myid] + nEqExtConstr + nEqBarConstr + nOfInstrConstr, (ONE / beta), knownTerms, 1); + + MPI_Barrier(MPI_COMM_WORLD); //only processor 0 accumulate and distribute the v array -#pragma omp for - for(i=0;i ZERO) + vVect[i] += vAuxVect[i]; + } + + nAstroElements = comlsqr.mapStar[myid][1] - comlsqr.mapStar[myid][0] + 1; + if (myid < nproc - 1) { - cblas_dscal ( nunkSplit, (ONE / alpha), vVect, 1 ); - cblas_dcopy ( nunkSplit, vVect, 1, wVect, 1 ); + nAstroElements = comlsqr.mapStar[myid][1] - comlsqr.mapStar[myid][0] + 1; + if (comlsqr.mapStar[myid][1] == comlsqr.mapStar[myid + 1][0]) + nAstroElements--; } + + double alphaLoc = 0; + if (nAstroPSolved) + alphaLoc = cblas_dnrm2(nAstroElements * comlsqr.nAstroPSolved, vVect, 1); + double alphaLoc2 = alphaLoc * alphaLoc; + if (myid == 0) + { + double alphaOther = cblas_dnrm2(comlsqr.nunkSplit - localAstroMax, &vVect[localAstroMax], 1); + double alphaOther2 = alphaOther * alphaOther; + alphaLoc2 = alphaLoc2 + alphaOther2; + } + double alphaGlob2 = 0; + MPI_Allreduce(&alphaLoc2, &alphaGlob2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + alpha = sqrt(alphaGlob2); +} + +if (alpha > ZERO) +{ + cblas_dscal(nunkSplit, (ONE / alpha), vVect, 1); + cblas_dcopy(nunkSplit, vVect, 1, wVect, 1); +} // printf("LSQR T4: PE=%d alpha=%15.12lf beta=%15.12lf arnorm=%15.12lf\n",myid,alpha,beta,arnorm); - arnorm = alpha * beta; - arnorm0 = arnorm; +arnorm = alpha * beta; +arnorm0 = arnorm; // printf("LSQR T5: PE=%d alpha=%15.12lf beta=%15.12lf arnorm=%15.12lf\n",myid,alpha,beta,arnorm); - if (arnorm == ZERO) goto goto_800; - rhobar = alpha; - phibar = beta; - bnorm = beta; - rnorm = beta; +if (arnorm == ZERO) + goto goto_800; +rhobar = alpha; +phibar = beta; +bnorm = beta; +rnorm = beta; - if (myid==0) +if (myid == 0) +{ + nout = fopen(lsqrOut, "a"); + if (nout != NULL) { - nout=fopen(lsqrOut,"a"); - if ( nout != NULL) { - if ( damped ) - fprintf(nout, fmt_1300); - else - fprintf(nout, fmt_1200); - } else { - printf("Error while opening %s (2)\n",lsqrOut); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - - test1 = ONE; - test2 = alpha / beta; - - if ( extra ) - fprintf(nout, fmt_1400); - - fprintf(nout, fmt_1550, itn, xSolution[0], rnorm, test1, test2); - fprintf(nout, fmt_1600); - fclose(nout); + if (damped) + fprintf(nout, fmt_1300); + else + fprintf(nout, fmt_1200); } + else + { + printf("Error while opening %s (2)\n", lsqrOut); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + + test1 = ONE; + test2 = alpha / beta; + + if (extra) + fprintf(nout, fmt_1400); + + fprintf(nout, fmt_1550, itn, xSolution[0], rnorm, test1, test2); + fprintf(nout, fmt_1600); + fclose(nout); +} // ================================================================== // CheckPointRestart // ================================================================== - - if(!noCPR) - if(myid==0) printf("PE=%d call restart setup\n",myid); - restartSetup( &itn, - knownTerms, - &beta, - &alpha, - vVect, - &anorm, - &rhobar, - &phibar, - wVect, - xSolution, - standardError, - &dnorm, - &sn2, - &cs2, - &z, - &xnorm1, - &res2, - &nstop, - comlsqr); - if(myid==0) printf("PE=%d end restart setup\n",myid); + +if (!noCPR) + if (myid == 0) + printf("PE=%d call restart setup\n", myid); +restartSetup(&itn, + knownTerms, + &beta, + &alpha, + vVect, + &anorm, + &rhobar, + &phibar, + wVect, + xSolution, + standardError, + &dnorm, + &sn2, + &cs2, + &z, + &xnorm1, + &res2, + &nstop, + comlsqr); +if (myid == 0) + printf("PE=%d end restart setup\n", myid); // ================================================================== // Main iteration loop. // ================================================================== - if(myid==0) - { - startCycleTime=time(NULL); - } -//////////////////////// START ITERATIONS - while (1) { - writeCPR=0; // No write CPR - cycleStartMpiTime=MPI_Wtime(); - itnTest=-1; - if(myid==0){ - fpItnLimitNew=NULL; - fpItnLimitNew=fopen("ITNLIMIT_NEW","r"); - if (fpItnLimitNew!=NULL) - { - fscanf(fpItnLimitNew, "%d\n",&itnTest); - if(itnTest>0) +if (myid == 0) +{ + startCycleTime = time(NULL); +} +//////////////////////// START ITERATIONS +while (1) +{ + writeCPR = 0; // No write CPR + cycleStartMpiTime = MPI_Wtime(); + itnTest = -1; + if (myid == 0) + { + fpItnLimitNew = NULL; + fpItnLimitNew = fopen("ITNLIMIT_NEW", "r"); + if (fpItnLimitNew != NULL) + { + fscanf(fpItnLimitNew, "%d\n", &itnTest); + if (itnTest > 0) { - comlsqr.itnLimit=itnTest; - if(itnlim != itnTest){ - itnlim=itnTest; - if(myid==0) printf("itnLimit forced with ITNLIMIT_NEW file to value %d\n",itnlim); + comlsqr.itnLimit = itnTest; + if (itnlim != itnTest) + { + itnlim = itnTest; + if (myid == 0) + printf("itnLimit forced with ITNLIMIT_NEW file to value %d\n", itnlim); } } fclose(fpItnLimitNew); - } } - - if(myid==0) + } + + if (myid == 0) + { + endCycleTime = time(NULL) - startCycleTime; + startCycleTime = time(NULL); + totTime = totTime + endCycleTime; + partialTime = partialTime + endCycleTime; + if (!noCPR) { - endCycleTime=time(NULL)-startCycleTime; - startCycleTime=time(NULL); - totTime=totTime+endCycleTime; - partialTime=partialTime+endCycleTime; - if(!noCPR){ - - if(partialTime>comlsqr.timeCPR*60) - { - writeCPR=1; - partialTime=0; - } - if(totTime+endCycleTime*2>comlsqr.timeLimit*60) writeCPR=2; - - if(CPRCount>0 && CPRCount%comlsqr.itnCPR==0) writeCPR=1; - if(CPRCount>0 && CPRCount==comlsqr.itnCPRstop) writeCPR=2; - if(CPRCount==itnlim) writeCPR=2; - fpCPRend=NULL; - fpCPRend=fopen("CPR_END","r"); - if (fpCPRend!=NULL) + + if (partialTime > comlsqr.timeCPR * 60) { - itnTest=-1; - fscanf(fpCPRend, "%d\n",&itnTest); - if(itnTest>0) + writeCPR = 1; + partialTime = 0; + } + if (totTime + endCycleTime * 2 > comlsqr.timeLimit * 60) + writeCPR = 2; + + if (CPRCount > 0 && CPRCount % comlsqr.itnCPR == 0) + writeCPR = 1; + if (CPRCount > 0 && CPRCount == comlsqr.itnCPRstop) + writeCPR = 2; + if (CPRCount == itnlim) + writeCPR = 2; + fpCPRend = NULL; + fpCPRend = fopen("CPR_END", "r"); + if (fpCPRend != NULL) + { + itnTest = -1; + fscanf(fpCPRend, "%d\n", &itnTest); + if (itnTest > 0) { - if(comlsqr.itnCPRend != itnTest){ - comlsqr.itnCPRend=itnTest; - printf("itnCPRend forced with CPR_END file to value %d\n",comlsqr.itnCPRend); + if (comlsqr.itnCPRend != itnTest) + { + comlsqr.itnCPRend = itnTest; + printf("itnCPRend forced with CPR_END file to value %d\n", comlsqr.itnCPRend); } } fclose(fpCPRend); } - if(comlsqr.itnCPRend>0 && comlsqr.itnCPRend<=itn) + if (comlsqr.itnCPRend > 0 && comlsqr.itnCPRend <= itn) { - writeCPR=2; - printf("itnCPRend condition triggered to value %d. writeCPR set to 2.\n",comlsqr.itnCPRend); + writeCPR = 2; + printf("itnCPRend condition triggered to value %d. writeCPR set to 2.\n", comlsqr.itnCPRend); } - }//if(!noCPR) - }//if(myid==0) - MPI_Barrier(MPI_COMM_WORLD); - MPI_Bcast( &writeCPR, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast( &itnlim, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast( &comlsqr.itnLimit, 1, MPI_INT, 0, MPI_COMM_WORLD); - if(myid==0) - { - printf("lsqr: Iteration number %d. Iteration seconds %ld. OmpSec %d Global Seconds %ld \n",itn,endCycleTime,ompSec,totTime); - if(writeCPR==1) printf("... writing CPR files\n"); - if(writeCPR==2) printf("... writing CPR files and return\n"); - } - - - if(writeCPR>0) - { - if(myid==0) - CPRtimeStart=time(NULL); - writeCheckPoint(itn, - knownTerms, - beta, - alpha, - vVect, - anorm, - rhobar, - phibar, - wVect, - xSolution, - standardError, - dnorm, - sn2, - cs2, - z, - xnorm1, - res2, - nstop, - comlsqr); - - if(myid==0){ - CPRtimeEnd=time(NULL)-CPRtimeStart; - totTime+=CPRtimeEnd; - partialTime+=CPRtimeEnd; - printf("CPR: itn=%d writing CPR seconds %ld. Global Seconds %ld \n",itn,CPRtimeEnd, totTime); - } - - if(writeCPR==2) - { - *istop_out = 1000; - *itn_out = itn; - if(nAstroPSolved) free(vAuxVect); - return; - } - if(myid==0) startCycleTime=time(NULL); - } + } //if(!noCPR) + } //if(myid==0) + MPI_Barrier(MPI_COMM_WORLD); + MPI_Bcast(&writeCPR, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&itnlim, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&comlsqr.itnLimit, 1, MPI_INT, 0, MPI_COMM_WORLD); + if (myid == 0) + { + printf("lsqr: Iteration number %d. Iteration seconds %ld. OmpSec %d Global Seconds %ld \n", itn, endCycleTime, ompSec, totTime); + if (writeCPR == 1) + printf("... writing CPR files\n"); + if (writeCPR == 2) + printf("... writing CPR files and return\n"); + } - itn = itn + 1; - comlsqr.itn=itn; - CPRCount++; -// ------------------------------------------------------------------ -// Perform the next step of the bidiagonalization to obtain the -// next beta, u, alpha, v. These satisfy the relations -// beta*u = A*v - alpha*u, -// alpha*v = A(transpose)*u - beta*v. -// ------------------------------------------------------------------ - cblas_dscal ( mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr, (- alpha), knownTerms, 1 ); - kAuxcopy=(double *) calloc(nEqExtConstr+nEqBarConstr+nOfInstrConstr, sizeof(double)); - for (int kk=0;kk 0) + { + if (myid == 0) + CPRtimeStart = time(NULL); + writeCheckPoint(itn, + knownTerms, + beta, + alpha, + vVect, + anorm, + rhobar, + phibar, + wVect, + xSolution, + standardError, + dnorm, + sn2, + cs2, + z, + xnorm1, + res2, + nstop, + comlsqr); + + if (myid == 0) + { + CPRtimeEnd = time(NULL) - CPRtimeStart; + totTime += CPRtimeEnd; + partialTime += CPRtimeEnd; + printf("CPR: itn=%d writing CPR seconds %ld. Global Seconds %ld \n", itn, CPRtimeEnd, totTime); } - aprod ( 1, m, n, vVect, knownTerms, systemMatrix, matrixIndex, instrCol,instrConstrIlung, comlsqr,&ompSec ); - kcopy=(double *) calloc(nEqExtConstr+nEqBarConstr+nOfInstrConstr, sizeof(double)); - MPI_Allreduce(&knownTerms[mapNoss[myid]],kcopy,nEqExtConstr+nEqBarConstr+nOfInstrConstr,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); - for(i=0;i ZERO) { - cblas_dscal ( mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr, (ONE / beta), knownTerms, 1 ); - cblas_dscal ( nunkSplit, (- beta), vVect, 1 ); -//only processor 0 accumulate and distribute the v array -#pragma omp for - for(i=0;i ZERO) { - vVect[localAstroMax+comlsqr.nAttParam+i]=dcopy[i]; - } - free(dcopy); - dcopy=(double *)calloc(comlsqr.nGlobalParam,sizeof(double)); - if(!dcopy) exit(err_malloc("dcopy",myid)); - mpi_allreduce(&vVect[localAstroMax+comlsqr.nAttParam+comlsqr.nInstrParam],dcopy,(long int) comlsqr.nGlobalParam,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); - for(i=0;i ZERO) { - cblas_dscal ( nunkSplit, (ONE / alpha), vVect, 1 ); - } + cblas_dscal(mapNoss[myid] + nEqExtConstr + nEqBarConstr + nOfInstrConstr, (ONE / beta), knownTerms, 1); + cblas_dscal(nunkSplit, (-beta), vVect, 1); +//only processor 0 accumulate and distribute the v array + // FV EDIT + //#pragma omp for + for (i = 0; i < localAstroMax; i++) + { + vAuxVect[i] = vVect[i]; + vVect[i] = 0; } - -// ------------------------------------------------------------------ -// Use a plane rotation to eliminate the damping parameter. -// This alters the diagonal (rhobar) of the lower-bidiagonal matrix. -// ------------------------------------------------------------------ - rhbar1 = rhobar; - if ( damped ) { - rhbar1 = d2norm( rhobar, damp ); - cs1 = rhobar / rhbar1; - sn1 = damp / rhbar1; - psi = sn1 * phibar; - phibar = cs1 * phibar; + if (myid != 0) + { + // FV EDIT + //#pragma omp for + for (i = localAstroMax; i < nunkSplit; i++) + vVect[i] = 0; } -// ------------------------------------------------------------------ -// Use a plane rotation to eliminate the subdiagonal element (beta) -// of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix. -// ------------------------------------------------------------------ - rho = d2norm( rhbar1, beta ); - cs = rhbar1 / rho; - sn = beta / rho; - theta = sn * alpha; - rhobar = - cs * alpha; - phi = cs * phibar; - phibar = sn * phibar; - tau = sn * phi; - -// ------------------------------------------------------------------ -// Update x, w and (perhaps) the standard error estimates. -// ------------------------------------------------------------------ - t1 = phi / rho; - t2 = - theta / rho; - t3 = ONE / rho; - dknorm = ZERO; - - for (i = 0; i < nAstroElements*comlsqr.nAstroPSolved; i++) { - t = wVect[i]; - t = (t3*t)*(t3*t); - dknorm = t + dknorm; + aprod(2, m, n, vVect, knownTerms, systemMatrix, matrixIndex, instrCol, instrConstrIlung, comlsqr, &ompSec ); + MPI_Barrier(MPI_COMM_WORLD); + //MPI_Finalize(); + //exit(0); + + //////////////// TEST vVect + /* + if(nproc==2) + { + printf("LSQR TP9 PE=%d itn=%d VrIdAstroPDim=%ld vVect[0-5]=%lf %lf %lf %lf %lf VrIdAstroPDim/2*5=%ld vVect[..+4]=%lf %lf %lf %lf %lf VrIdAstroPDim-1*5=%ld vVect[...+5]=%lf %lf %lf %lf %lf\n",myid,itn,VrIdAstroPDim,vVect[0],vVect[1],vVect[2],vVect[3],vVect[4],(VrIdAstroPDim/2)*5,vVect[(VrIdAstroPDim/2)*5],vVect[(VrIdAstroPDim/2)*5+1],vVect[(VrIdAstroPDim/2)*5+2],vVect[(VrIdAstroPDim/2)*5+3],vVect[(VrIdAstroPDim/2)*5+4],(VrIdAstroPDim-1)*5,vVect[(VrIdAstroPDim-1)*5],vVect[(VrIdAstroPDim-1)*5+1],vVect[(VrIdAstroPDim-1)*5+2],vVect[(VrIdAstroPDim-1)*5+3],vVect[(VrIdAstroPDim-1)*5+4]); } - dknormSum=0; - MPI_Allreduce(&dknorm,&dknormSum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); - dknorm=dknormSum; - - if ( wantse ) { -#pragma omp for - for (i = 0; i < localAstro; i++) { - t = wVect[i]; - xSolution[i] = t1*t + xSolution[i]; - wVect[i] = t2*t + vVect[i]; - t = (t3*t)*(t3*t); - standardError[i] = t + standardError[i]; - } -#pragma omp for - for (i = localAstroMax; i < localAstroMax+other; i++) { - t = wVect[i]; - xSolution[i] = t1*t + xSolution[i]; - wVect[i] = t2*t + vVect[i]; - t = (t3*t)*(t3*t); - standardError[i] = t + standardError[i]; - } - for (i = localAstroMax; i < localAstroMax+other; i++) { - t = wVect[i]; - t = (t3*t)*(t3*t); - dknorm = t + dknorm; - } + if(nproc==1) + { + printf("LSQR TP9 PE=%d itn=%d VrIdAstroPDim=%ld vVect[0-5]=%lf %lf %lf %lf %lf VrIdAstroPDim/2*5=%ld vVect[..+4]=%lf %lf %lf %lf %lf VrIdAstroPDim-1*5=%ld vVect[...+5]=%lf %lf %lf %lf %lf\n",myid,itn,VrIdAstroPDim,vVect[0],vVect[1],vVect[2],vVect[3],vVect[4],(VrIdAstroPDim/2)*5,vVect[(VrIdAstroPDim/2)*5],vVect[(VrIdAstroPDim/2)*5+1],vVect[(VrIdAstroPDim/2)*5+2],vVect[(VrIdAstroPDim/2)*5+3],vVect[(VrIdAstroPDim/2)*5+4],(VrIdAstroPDim-1)*5,vVect[(VrIdAstroPDim-1)*5],vVect[(VrIdAstroPDim-1)*5+1],vVect[(VrIdAstroPDim-1)*5+2],vVect[(VrIdAstroPDim-1)*5+3],vVect[(VrIdAstroPDim-1)*5+4]); } - else { -#pragma omp for - for (i = 0; i < localAstro; i++) { - t = wVect[i]; - xSolution[i] = t1*t + xSolution[i]; - wVect[i] = t2*t + vVect[i]; - } -#pragma omp for - for (i = localAstroMax; i < localAstroMax+other; i++) { - t = wVect[i]; - xSolution[i] = t1*t + xSolution[i]; - wVect[i] = t2*t + vVect[i]; - } - for (i = localAstroMax; i < localAstroMax+other; i++) { - t = wVect[i]; - dknorm = (t3*t)*(t3*t) + dknorm; - } + */ + ////////////////////////// + + double *dcopy; + dcopy = (double *)calloc(comlsqr.nAttParam, sizeof(double)); + if (!dcopy) + exit(err_malloc("dcopy", myid)); + mpi_allreduce(&vVect[localAstroMax], dcopy, (long int)comlsqr.nAttParam, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + // FV EDIT + //#pragma omp for + for (i = 0; i < comlsqr.nAttParam; i++) + { + vVect[localAstroMax + i] = dcopy[i]; + } + free(dcopy); + dcopy = (double *)calloc(comlsqr.nInstrParam, sizeof(double)); + if (!dcopy) + exit(err_malloc("dcopy", myid)); + mpi_allreduce(&vVect[localAstroMax + comlsqr.nAttParam], dcopy, (long int)comlsqr.nInstrParam, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + // FV EDIT + //#pragma omp for + for (i = 0; i < comlsqr.nInstrParam; i++) + { + vVect[localAstroMax + comlsqr.nAttParam + i] = dcopy[i]; + } + free(dcopy); + dcopy = (double *)calloc(comlsqr.nGlobalParam, sizeof(double)); + if (!dcopy) + exit(err_malloc("dcopy", myid)); + mpi_allreduce(&vVect[localAstroMax + comlsqr.nAttParam + comlsqr.nInstrParam], dcopy, (long int)comlsqr.nGlobalParam, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + for (i = 0; i < comlsqr.nGlobalParam; i++) + { + vVect[localAstroMax + comlsqr.nAttParam + comlsqr.nInstrParam + i] = dcopy[i]; + } + free(dcopy); + if (nAstroPSolved) + SumCirc(vVect, comlsqr); + // FV EDIT + //#pragma omp for + for (i = 0; i < localAstroMax; i++) + { + vVect[i] += vAuxVect[i]; } -// ------------------------------------------------------------------ -// Monitor the norm of d_k, the update to x. -// dknorm = norm( d_k ) -// dnorm = norm( D_k ), where D_k = (d_1, d_2, ..., d_k ) -// dxk = norm( phi_k d_k ), where new x = x_k + phi_k d_k. -// ------------------------------------------------------------------ - dknorm = sqrt( dknorm ); - dnorm = d2norm( dnorm, dknorm ); - dxk = fabs( phi * dknorm ); - if (dxmax < dxk ) { - dxmax = dxk; - maxdx = itn; + double alphaLoc = 0; + if (nAstroPSolved) + alphaLoc = cblas_dnrm2(nAstroElements * comlsqr.nAstroPSolved, vVect, 1); + double alphaLoc2 = alphaLoc * alphaLoc; + if (myid == 0) + { + double alphaOther = cblas_dnrm2(comlsqr.nunkSplit - localAstroMax, &vVect[localAstroMax], 1); + double alphaOther2 = alphaOther * alphaOther; + alphaLoc2 = alphaLoc2 + alphaOther2; } + double alphaGlob2 = 0; + MPI_Allreduce(&alphaLoc2, &alphaGlob2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); -// ------------------------------------------------------------------ -// Use a plane rotation on the right to eliminate the -// super-diagonal element (theta) of the upper-bidiagonal matrix. -// Then use the result to estimate norm(x). -// ------------------------------------------------------------------ - delta = sn2 * rho; - gambar = - cs2 * rho; - rhs = phi - delta * z; - zbar = rhs / gambar; - xnorm = d2norm( xnorm1, zbar ); - gamma = d2norm( gambar, theta ); - cs2 = gambar / gamma; - sn2 = theta / gamma; - z = rhs / gamma; - xnorm1 = d2norm( xnorm1, z ); - -// ------------------------------------------------------------------ -// Test for convergence. -// First, estimate the norm and condition of the matrix Abar, -// and the norms of rbar and Abar(transpose)*rbar. -// ------------------------------------------------------------------ - acond = anorm * dnorm; - res2 = d2norm( res2 , psi ); - rnorm = d2norm( res2 , phibar ); - arnorm = alpha * fabs( tau ); - -// Now use these norms to estimate certain other quantities, -// some of which will be small near a solution. - - alfopt = sqrt( rnorm / (dnorm * xnorm) ); - test1 = rnorm / bnorm; - test2 = ZERO; - if (rnorm > ZERO) test2 = arnorm / (anorm * rnorm); -// if (arnorm0 > ZERO) test2 = arnorm / arnorm0; //(Michael Friedlander's modification) - test3 = ONE / acond; - t1 = test1 / (ONE + anorm * xnorm / bnorm); - rtol = btol + atol * anorm * xnorm / bnorm; - -// The following tests guard against extremely small values of -// atol, btol or ctol. (The user may have set any or all of -// the parameters atol, btol, conlim to zero.) -// The effect is equivalent to the normal tests using -// atol = relpr, btol = relpr, conlim = 1/relpr. - - t3 = ONE + test3; - t2 = ONE + test2; - t1 = ONE + t1; -// printf("LSQR TP8 PE=%d itn=%d test1=%18.16lf, rnorm=%lf bnorm=%lf anorm=%lf xnorm=%lf rtol=%lf t1=%18.16lf\n",myid,itn,test1,rnorm,bnorm,anorm,xnorm,rtol,t1); - - if (itn >= itnlim) istop = 5; - if (t3 <= ONE ) istop = 4; - if (t2 <= ONE ) istop = 2; - if (t1 <= ONE ) istop = 1; -// if (t1 <= ONE ) printf("PE=%d t1=%lf\n",myid,t1); - -// Allow for tolerances set by the user. - - if (test3 <= ctol) istop = 4; - if (test2 <= atol) istop = 2; - if (test1 <= rtol) istop = 1; //(Michael Friedlander had this commented out) -// if (test1 <= rtol) printf("PE=%d test1=%lf\n",myid,test1); - -// ------------------------------------------------------------------ -// See if it is time to print something. -// ------------------------------------------------------------------ -// if (nout == NULL ) goto goto_600; // Delete for buffer flush modification??? TBV - if (n <= 40 ) goto goto_400; - if (itn <= 10 ) goto goto_400; - if (itn >= itnlim-10) goto goto_400; - if (itn % 10 == 0 ) goto goto_400; - if (test3 <= 2.0*ctol) goto goto_400; - if (test2 <= 10.0*atol) goto goto_400; - if (test1 <= 10.0*rtol) goto goto_400; - if (istop != 0 ) goto goto_400; - goto goto_600; - -// Print a line for this iteration. -// "extra" is for experimental purposes. - - goto_400: - if(!comlsqr.Test) - if(myid==0) { - nout=fopen(lsqrOut,"a"); - if ( nout != NULL) { - if ( extra ) { - fprintf(nout, fmt_1500_extra, - itn, xSolution[0], rnorm, test1, test2, anorm, - acond, phi, dknorm, dxk, alfopt); - } else { - fprintf(nout, fmt_1500, - itn, xSolution[0], rnorm, test1, test2, anorm, acond); - } - if (itn % 10 == 0) fprintf(nout, fmt_1600); - } else { - printf("Error while opening %s (3)\n",lsqrOut); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - fclose(nout); + alpha = sqrt(alphaGlob2); + if (alpha > ZERO) + { + cblas_dscal(nunkSplit, (ONE / alpha), vVect, 1); } -// ------------------------------------------------------------------ -// Stop if appropriate. -// The convergence criteria are required to be met on nconv -// consecutive iterations, where nconv is set below. -// Suggested value: nconv = 1, 2 or 3. -// ------------------------------------------------------------------ - goto_600: - if (istop == 0) { - nstop = 0; + } + + // ------------------------------------------------------------------ + // Use a plane rotation to eliminate the damping parameter. + // This alters the diagonal (rhobar) of the lower-bidiagonal matrix. + // ------------------------------------------------------------------ + rhbar1 = rhobar; + if (damped) + { + rhbar1 = d2norm(rhobar, damp); + cs1 = rhobar / rhbar1; + sn1 = damp / rhbar1; + psi = sn1 * phibar; + phibar = cs1 * phibar; + } + + // ------------------------------------------------------------------ + // Use a plane rotation to eliminate the subdiagonal element (beta) + // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix. + // ------------------------------------------------------------------ + rho = d2norm(rhbar1, beta); + cs = rhbar1 / rho; + sn = beta / rho; + theta = sn * alpha; + rhobar = -cs * alpha; + phi = cs * phibar; + phibar = sn * phibar; + tau = sn * phi; + + // ------------------------------------------------------------------ + // Update x, w and (perhaps) the standard error estimates. + // ------------------------------------------------------------------ + t1 = phi / rho; + t2 = -theta / rho; + t3 = ONE / rho; + dknorm = ZERO; + + for (i = 0; i < nAstroElements * comlsqr.nAstroPSolved; i++) + { + t = wVect[i]; + t = (t3 * t) * (t3 * t); + dknorm = t + dknorm; + } + dknormSum = 0; + MPI_Allreduce(&dknorm, &dknormSum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + dknorm = dknormSum; + + if (wantse) + { + // FV EDIT + //#pragma omp for + for (i = 0; i < localAstro; i++) + { + t = wVect[i]; + xSolution[i] = t1 * t + xSolution[i]; + wVect[i] = t2 * t + vVect[i]; + t = (t3 * t) * (t3 * t); + standardError[i] = t + standardError[i]; } - else { - nconv = 1; - nstop = nstop + 1; - if (nstop < nconv && itn < itnlim) istop = 0; + // FV EDIT + //#pragma omp for + for (i = localAstroMax; i < localAstroMax + other; i++) + { + t = wVect[i]; + xSolution[i] = t1 * t + xSolution[i]; + wVect[i] = t2 * t + vVect[i]; + t = (t3 * t) * (t3 * t); + standardError[i] = t + standardError[i]; } - - if(comlsqr.Test) + for (i = localAstroMax; i < localAstroMax + other; i++) { - if(itn ZERO) + test2 = arnorm / (anorm * rnorm); + // if (arnorm0 > ZERO) test2 = arnorm / arnorm0; //(Michael Friedlander's modification) + test3 = ONE / acond; + t1 = test1 / (ONE + anorm * xnorm / bnorm); + rtol = btol + atol * anorm * xnorm / bnorm; + + // The following tests guard against extremely small values of + // atol, btol or ctol. (The user may have set any or all of + // the parameters atol, btol, conlim to zero.) + // The effect is equivalent to the normal tests using + // atol = relpr, btol = relpr, conlim = 1/relpr. + + t3 = ONE + test3; + t2 = ONE + test2; + t1 = ONE + t1; + // printf("LSQR TP8 PE=%d itn=%d test1=%18.16lf, rnorm=%lf bnorm=%lf anorm=%lf xnorm=%lf rtol=%lf t1=%18.16lf\n",myid,itn,test1,rnorm,bnorm,anorm,xnorm,rtol,t1); + + if (itn >= itnlim) + istop = 5; + if (t3 <= ONE) + istop = 4; + if (t2 <= ONE) + istop = 2; + if (t1 <= ONE) + istop = 1; + // if (t1 <= ONE ) printf("PE=%d t1=%lf\n",myid,t1); + + // Allow for tolerances set by the user. + + if (test3 <= ctol) + istop = 4; + if (test2 <= atol) + istop = 2; + if (test1 <= rtol) + istop = 1; //(Michael Friedlander had this commented out) + // if (test1 <= rtol) printf("PE=%d test1=%lf\n",myid,test1); + + // ------------------------------------------------------------------ + // See if it is time to print something. + // ------------------------------------------------------------------ + // if (nout == NULL ) goto goto_600; // Delete for buffer flush modification??? TBV + if (n <= 40) + goto goto_400; + if (itn <= 10) + goto goto_400; + if (itn >= itnlim - 10) + goto goto_400; + if (itn % 10 == 0) + goto goto_400; + if (test3 <= 2.0 * ctol) + goto goto_400; + if (test2 <= 10.0 * atol) + goto goto_400; + if (test1 <= 10.0 * rtol) + goto goto_400; + if (istop != 0) + goto goto_400; + goto goto_600; + + // Print a line for this iteration. + // "extra" is for experimental purposes. + +goto_400: + if (!comlsqr.Test) + if (myid == 0) + { + nout = fopen(lsqrOut, "a"); + if (nout != NULL) + { + if (extra) + { + fprintf(nout, fmt_1500_extra, + itn, xSolution[0], rnorm, test1, test2, anorm, + acond, phi, dknorm, dxk, alfopt); + } + else + { + fprintf(nout, fmt_1500, + itn, xSolution[0], rnorm, test1, test2, anorm, acond); + } + if (itn % 10 == 0) + fprintf(nout, fmt_1600); + } + else + { + printf("Error while opening %s (3)\n", lsqrOut); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + fclose(nout); + } + // ------------------------------------------------------------------ + // Stop if appropriate. + // The convergence criteria are required to be met on nconv + // consecutive iterations, where nconv is set below. + // Suggested value: nconv = 1, 2 or 3. + // ------------------------------------------------------------------ +goto_600: + if (istop == 0) + { + nstop = 0; + } + else + { + nconv = 1; + nstop = nstop + 1; + if (nstop < nconv && itn < itnlim) + istop = 0; + } + + if (comlsqr.Test) + { + if (itn < comlsqr.itnLimit) + istop = 0; + else + (istop = 5); + } + cycleEndMpiTime = MPI_Wtime(); + // printf("lsqr: PE=%d MPI-iteration number %d. Iteration seconds %f\n",myid,itn,cycleEndMpiTime-cycleStartMpiTime); + + if (istop != 0) + break; +} // ================================================================== // End of iteration loop. // ================================================================== // Finish off the standard error estimates. - if ( wantse ) { - t = ONE; - if (m > n) t = m - n; - if ( damped ) t = m; - t = rnorm / sqrt( t ); - - for (i = 0; i < nunkSplit; i++) - standardError[i] = t * sqrt( standardError[i] ); - - } +if (wantse) +{ + t = ONE; + if (m > n) + t = m - n; + if (damped) + t = m; + t = rnorm / sqrt(t); + + for (i = 0; i < nunkSplit; i++) + standardError[i] = t * sqrt(standardError[i]); +} // Decide if istop = 2 or 3. // Print the stopping condition. - goto_800: - if (damped && istop == 2) istop = 3; - if (myid == 0) { - nout=fopen(lsqrOut,"a"); - if (nout != NULL) { - fprintf(nout, fmt_2000, - exitLsqr, istop, itn, - exitLsqr, anorm, acond, - exitLsqr, bnorm, xnorm, - exitLsqr, rnorm, arnorm); - fprintf(nout, fmt_2100, - exitLsqr, dxmax, maxdx, - exitLsqr, dxmax/(xnorm + 1.0e-20)); - fprintf(nout, fmt_3000, - exitLsqr, msg[istop]); - } else { - printf("Error while opening %s (4)\n",lsqrOut); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - fclose(nout); - } +goto_800 : if (damped && istop == 2) istop = 3; +if (myid == 0) +{ + nout = fopen(lsqrOut, "a"); + if (nout != NULL) + { + fprintf(nout, fmt_2000, + exitLsqr, istop, itn, + exitLsqr, anorm, acond, + exitLsqr, bnorm, xnorm, + exitLsqr, rnorm, arnorm); + fprintf(nout, fmt_2100, + exitLsqr, dxmax, maxdx, + exitLsqr, dxmax / (xnorm + 1.0e-20)); + fprintf(nout, fmt_3000, + exitLsqr, msg[istop]); + } + else + { + printf("Error while opening %s (4)\n", lsqrOut); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + fclose(nout); +} // Assign output variables from local copies. - *istop_out = istop; - *itn_out = itn; - *anorm_out = anorm; - *acond_out = acond; - *rnorm_out = rnorm; - *arnorm_out = test2; - *xnorm_out = xnorm; - - if(nAstroPSolved) free(vAuxVect); - return; +*istop_out = istop; +*itn_out = itn; +*anorm_out = anorm; +*acond_out = acond; +*rnorm_out = rnorm; +*arnorm_out = test2; +*xnorm_out = xnorm; + +if (nAstroPSolved) + free(vAuxVect); +return; } - - diff --git a/lsqrblas.c b/lsqrblas.c index 7025ed1ecd1e780d50490aa6495cf9d6a6b49f7d..1f3f715ab9011889c768f11e8deee4debbd39512 100644 --- a/lsqrblas.c +++ b/lsqrblas.c @@ -26,6 +26,10 @@ */ #include #include +#define MIN(a,b) (((a)<(b))?(a):(b)) +#define MAX(a,b) (((a)>(b))?(a):(b)) + + /*! \file @@ -95,7 +99,7 @@ cblas_daxpy( const int N, const double alpha, const double *X, */ void cblas_dcopy( const long int N, const double *X, - const int incX, double *Y, const int incY) + const int incX, double *Y, const int incY, int ntasks) { long int i,ix,iy; // int ix = OFFSET(N, incX); @@ -107,17 +111,24 @@ cblas_dcopy( const long int N, const double *X, else iy=((N) - 1) * (-(incY)); -#pragma omp parallel shared(ix, iy) +//#pragma omp parallel shared(ix, iy) +//#pragma omp task shared(ix, iy) { + + long ixThread,iyThread; long ixThreadOffset=ix, iyThreadOffset=iy; -#pragma omp for +//#pragma omp for for (i = 0; i < N; i++) { + //printf("Inside task %d\n", i); ixThread=ixThreadOffset+incX*i; iyThread=iyThreadOffset+incY*i; Y[iyThread] = X[ixThread]; } + + }//end pragma +//#pragma omp taskwait } /*! @@ -157,7 +168,7 @@ cblas_ddot( const int N, const double *X, */ double -cblas_dnrm2( const long int N, const double *X, const int incX) +cblas_dnrm2( const long int N, const double *X, const int incX, int ntasks) { double scale = 0.0, @@ -190,9 +201,9 @@ cblas_dnrm2( const long int N, const double *X, const int incX) return scale * sqrt(ssq); } - +//FV_ not used ?? double -cblas_dnrm21( const long int N, const double *X, const int incX) +cblas_dnrm21( const long int N, const double *X, const int incX, int ntasks=0) { double ssq[100], ssqFinal=0; @@ -204,21 +215,25 @@ cblas_dnrm21( const long int N, const double *X, const int incX) if (N <= 0 || incX <= 0) return 0; else if (N == 1) return fabs(X[0]); -#pragma omp parallel private(tid,nthreads,ix) + +////#pragma omp parallel private(tid,nthreads,ix) { + + /* #ifdef OMP - tid = omp_get_thread_num(); + tid = omp_get_thread_num(); nthreads = omp_get_num_threads(); if(nthreads>100) exit(1); #endif - + */ -#pragma omp for +////#pragma omp for for (long i = 0; i < N; i++) { ix=incX*i; ssq[tid] += X[ix]*X[ix]; } -} // pragma +} +// pragma for(int j=0;j<100;j++) ssqFinal+=ssq[j]; @@ -226,7 +241,7 @@ for(int j=0;j<100;j++) return sqrt(ssqFinal); } - +//FV_ not used ?? double cblas_dnrm22( const long int N, const double *X, const int incX) { @@ -243,14 +258,14 @@ cblas_dnrm22( const long int N, const double *X, const int incX) if (N <= 0 || incX <= 0) return 0; else if (N == 1) return fabs(X[0]); -#pragma omp parallel private(ix) -{ -#ifdef OMP - tid = omp_get_thread_num(); - nthreads = omp_get_num_threads(); - if(nthreads>100) exit(1); -#endif -#pragma omp for +///#pragma omp parallel private(ix) +////{ +////#ifdef OMP +//// tid = omp_get_thread_num(); +//// nthreads = omp_get_num_threads(); +//// if(nthreads>100) exit(1); +////#endif +////#pragma omp for for (i = 0; i < N; i++) { ix=incX*i; const double x = X[ix]; @@ -269,7 +284,7 @@ cblas_dnrm22( const long int N, const double *X, const int incX) } resultSqr[tid]=(scale[tid] * sqrt(ssq[tid]))*(scale[tid] * sqrt(ssq[tid])); -}// pragma +////}// pragma double resultFinal=0.0; for(int j=0;j<100;j++) resultFinal+=resultSqr[j]; @@ -301,6 +316,7 @@ void cblas_dscal(const long int N, const double alpha, double *X, const int incX } } +//FV_ not used ?? void cblas_dscal1(const long int N, const double alpha, double *X, const int incX) { @@ -312,14 +328,14 @@ cblas_dscal1(const long int N, const double alpha, double *X, const int incX) if (incX > 0) ix=0; else ix=((N) - 1) * (-(incX)); -#pragma omp parallel shared(ix) -{ +////#pragma omp parallel shared(ix) +/////{ long ixThread; long ixThreadOffset=ix; -#pragma omp for +/////#pragma omp for for (i = 0; i < N; i++) { ixThread=ixThreadOffset+incX*i; X[ixThread] *= alpha; } -}// end pragma +/////}// end pragma } diff --git a/memRequest.c b/memRequest.c old mode 100755 new mode 100644 diff --git a/solvergaiaSim.c b/solvergaiaSim.c index d001f9ad8618aa8f4c38c26d998eaa5d63086274..02dcb17100d311d03a8b7041a9687817eeaf5887 100644 --- a/solvergaiaSim.c +++ b/solvergaiaSim.c @@ -5,9 +5,9 @@ M.G. Lattanzi, A. Vecchiato, B. Bucciarelli (23 May 1996) A. Vecchiato, R. Morbidelli for the Fortran 90 version with dynamical memory allocation (21 March 2005) - + A. Vecchiato, June 16, 2008 - + Version 2.1 , Feb 8, 2012 Version history: - version 2.1, Feb 8 2012 - Added debug mode @@ -16,8 +16,6 @@ - version 5.0 - May 2013 from Ugo Becciani and Alberto Vecchiato */ - - //#include #include #include @@ -34,7 +32,6 @@ #define MAX_CONSTR 1000 - long instr_hash(int FoV, int CCD, int PixelColumn, int TimeInterval); long randlong(long max); long randlong1(long min, long max); @@ -42,505 +39,581 @@ int randint(int max); int randint1(int min, int max); int fill_extract(long *values, long *pos_min, long pos_max, long *number); - - /* Start of main program */ -int main(int argc, char **argv) { - int debugMode,wrsol; +int main(int argc, char **argv) +{ + int debugMode, wrsol; int i; - - int idtest=0, precond=1; + + int idtest = 0, precond = 1; double srIDtest, pert; - - int inputDirOpt=0, outputDirOpt=0, inputDirLen; - char inputDir[1024]="", outputDir[1024]="",wrfileDir[1024]=""; + + int inputDirOpt = 0, outputDirOpt = 0, inputDirLen; + char inputDir[1024] = "", outputDir[1024] = "", wrfileDir[1024] = ""; char wpath[1024]; - size_t sizePath=1020; - + size_t sizePath = 1020; + char filenameSolProps[150]; char filenameSolPropsFinal[150]; - char filenameAstroResults[150]; /* file storing the Astrometric Parameters */ - char filenameAttResults[150]; /* file storing the Attitude Parameters */ - char filenameInstrResults[150]; /* file storing the Instrument Parameters */ + char filenameAstroResults[150]; /* file storing the Astrometric Parameters */ + char filenameAttResults[150]; /* file storing the Attitude Parameters */ + char filenameInstrResults[150]; /* file storing the Instrument Parameters */ char filenameGlobalResults[150]; /* file storing the Global Parameters */ - + long ii, jj, kk; long sphereId; // Id of the sphere to be solved - long idum; // variable to initialize the random number generator aggiunta la variabile per inizializzare ran2() - + long idum; // variable to initialize the random number generator aggiunta la variabile per inizializzare ran2() + // LSQR input parameters - long itnlim; // maximum number of iteration allowed for the LSQR convergence + long itnlim; // maximum number of iteration allowed for the LSQR convergence double damp, atol, btol, conlim; // other LSQR input parameters - double aTol; //read by command line, overrides atol if >=0 - + double aTol; //read by command line, overrides atol if >=0 + // LSQR output parameters - int istop; // LSQR stopping condition - int itn; // LSQR iteration number + int istop; // LSQR stopping condition + int itn; // LSQR iteration number double anorm, acond, rnorm, arnorm, xnorm; // other LSQR output parameters - + // Properties of the system to be solved // Astrometric parameters: - long nStar; // number of stars - short nAstroP; // number of astrometric parameters for each observation - - short nAstroPSolved,nInstrPSolved; // number of astrometric and instrumental parameters taken as unknowns in the system of equations - int *mapAstroP; // indeces of the solved astrometric parameters - int lsInstrFlag,ssInstrFlag,nuInstrFlag,maInstrFlag; - int nElemIC=0,nOfInstrConstr=0; - int nOfInstrConstrLSAL=0, nElemICLSAL=0, nOfInstrConstrLSAC=0, nElemICLSAC=0 , nOfInstrConstrSS=0 ,nElemICSS=0; + long nStar; // number of stars + short nAstroP; // number of astrometric parameters for each observation + + short nAstroPSolved, nInstrPSolved; // number of astrometric and instrumental parameters taken as unknowns in the system of equations + int *mapAstroP; // indeces of the solved astrometric parameters + int lsInstrFlag, ssInstrFlag, nuInstrFlag, maInstrFlag; + int nElemIC = 0, nOfInstrConstr = 0; + int nOfInstrConstrLSAL = 0, nElemICLSAL = 0, nOfInstrConstrLSAC = 0, nElemICLSAC = 0, nOfInstrConstrSS = 0, nElemICSS = 0; // Attitude parameters: - - - long nDegFreedomAtt=0; // number of degrees of freedom for each attitude axis - long cCDLSAACZP=14; // CCD light sensitive area AC zero point - short nAttAxes=0, nAttParAxis; // number of attitude axes to be reconstructed, number of non-zero attitude coefficients per obs. & axis - short nAttP=0; // number of attitude parameters for each observation + + long nDegFreedomAtt = 0; // number of degrees of freedom for each attitude axis + long cCDLSAACZP = 14; // CCD light sensitive area AC zero point + short nAttAxes = 0, nAttParAxis; // number of attitude axes to be reconstructed, number of non-zero attitude coefficients per obs. & axis + short nAttP = 0; // number of attitude parameters for each observation // Instrument parameters: long instrSetUp; // number coding the set up of the instrument in terms of: # of FoVs, # of CCDs, # of pixel columns, # of time intervals short nInstrP; // number of instrument parameters for each observation // Global parameters: - short nGlobP; // number of global parameters - - long nobs; // number of observations - long nConstrLong, nConstrLat, nConstrMuLong, nConstrMuLat; // number of constraints for the astrometric parameters + short nGlobP; // number of global parameters + + long nobs; // number of observations + long nConstrLong, nConstrLat, nConstrMuLong, nConstrMuLat; // number of constraints for the astrometric parameters long *constrLongId, *constrLatId, *constrMuLongId, *constrMuLatId; // arrays with the gsrIDs of the constrained sources - double *constrLongW, *constrLatW, *constrMuLongW, *constrMuLatW; // arrays with the weights of the constraints - double extConstrW; // weight of the null space constraint equations - double barConstrW; // weight of the baricentric constraint equations - const double attExtConstrFact=-0.5; // the coefficients of the attitude part of the null space constraint equations must be multiplied by a factor -1/2 - + double *constrLongW, *constrLatW, *constrMuLongW, *constrMuLatW; // arrays with the weights of the constraints + double extConstrW; // weight of the null space constraint equations + double barConstrW; // weight of the baricentric constraint equations + const double attExtConstrFact = -0.5; // the coefficients of the attitude part of the null space constraint equations must be multiplied by a factor -1/2 + // Variables to be read from the GsrSystemRows - + // Internal variables long nAstroParam; - int nAttParam, nInstrParam, nGlobalParam,nInstrParamTot; // total number of unknowns for the Astrometric, Attitude, Instrument, and Global parameters respectively - int nparam; // total number of unknowns + int nAttParam, nInstrParam, nGlobalParam, nInstrParamTot; // total number of unknowns for the Astrometric, Attitude, Instrument, and Global parameters respectively + int nparam; // total number of unknowns long nunk, nunkSplit, nConstr; - long ncoeff, ielem, global_ielem, precnofrows,global_precnofrows, ncolumn, nrowsToRead; + long ncoeff, ielem, global_ielem, precnofrows, global_precnofrows, ncolumn, nrowsToRead; long offsetAttParam, offsetInstrParam, offsetGlobParam, VroffsetAttParam; // offests for the Attitude, Instrument, and Global columns of the coefficient matrix - unsigned long int totmem, nElements; // total required memory and temporary variable to store the memory being allocated - long VrIdAstroPDimMax=0; // - long VrIdAstroPDim=0; // + unsigned long int totmem, nElements; // total required memory and temporary variable to store the memory being allocated + long VrIdAstroPDimMax = 0; // + long VrIdAstroPDim = 0; // long VrIdAstroPDimRecv; int offset; long nObsxStar, nobsOri; - int nfileProc=3; - int addElementAtt=0; - int addElementextStar=0; - int addElementbarStar=0; - int nOfElextObs=0; - int nOfElBarObs=0; + int nfileProc = 3; + int addElementAtt = 0; + int addElementextStar = 0; + int addElementbarStar = 0; + int nOfElextObs = 0; + int nOfElBarObs = 0; double *attNS; - // Array arguments of the LSQR function - double *knownTerms, *vVect, *wVect, *xSolution, *standardError; // arrays containing the known terms (knownterms), the bidiagonalization (wVect, wVect), the unknowns (xSolution), and the estimated variances (standardError) - long int *matrixIndex; // ia[i] contains the column number of the coefficient systemMatrix[i] + double *knownTerms, *vVect, *wVect, *xSolution, *standardError; + + double **vVect_aux_AttP; + double **vVect_aux_InstP; + + // arrays containing the known terms (knownterms), the bidiagonalization (wVect, wVect), the unknowns (xSolution), and the estimated variances (standardError) + long int *matrixIndex; // ia[i] contains the column number of the coefficient systemMatrix[i] int *instrConst; - int *instrCol; // columns on each observation for instumental Parameters - int * instrConstrIlung; // vector containing the length of each instrConstr eq. + int *instrCol; // columns on each observation for instumental Parameters + int *instrConstrIlung; // vector containing the length of each instrConstr eq. double *systemMatrix, *preCondVect; // array of the non-zero coefficients of the system and of the column normalization respectively - int wgInstrCoeff=1; - + int wgInstrCoeff = 1; + ///////////////////// // Arrays containing the solution coming from the LSQR - double *xAstro, *standardErrorAstro; // solution and standard errors for the Astrometric parameters - double *xAtt, *standardErrorAtt; // solution and standard errors for the Attitude parameters - double *xInstr, *standardErrorInstr; // solution and standard errors for the Instrument parameters + double *xAstro, *standardErrorAstro; // solution and standard errors for the Astrometric parameters + double *xAtt, *standardErrorAtt; // solution and standard errors for the Attitude parameters + double *xInstr, *standardErrorInstr; // solution and standard errors for the Instrument parameters double *xGlobal, *standardErrorGlobal; // solution and standard errors for the Global parameters - - + // file pointers - FILE *fpSolProps,*fpSolPropsFinal,*fpSolPropsFileBin; + FILE *fpSolProps, *fpSolPropsFinal, *fpSolPropsFileBin; FILE *fpMI, *fpSM, *fpII, *fpKT, *fpCPR; FILE *fpAstroR, *fpAttR, *fpInstrR, *fpGlobR, *ffcont; - + time_t seconds[10], tot_sec[10]; seconds[0] = time(NULL); double timeToReadFiles; - + //MPI definitions - int nproc, myid, namelen; + int nproc, myid, namelen; char processor_name[MPI_MAX_PROCESSOR_NAME]; - double startTime,endTime; - + double startTime, endTime; + // Distributed array maps - long int * mapNcoeff, * mapNoss; - long int mapNcoeffBefore, mapNossBefore; - long int mapNcoeffAfter, mapNossAfter; + long int *mapNcoeff, *mapNoss; + long int mapNcoeffBefore, mapNossBefore; + long int mapNcoeffAfter, mapNossAfter; // struct to communicate with lsqr - int timeCPR, timeLimit,itnCPR,itnCPRstop=1000000; + int timeCPR, timeLimit, itnCPR, itnCPRstop = 1000000; int itnLimit; //read by command line, overrides itnlim if >0 struct comData comlsqr; // serach for map - long lastStar=-1; - long firstStar=-1; - long lastStarConstr=-1; - long firstStarConstr=-1; - int seqStar=1; - int withFile=1; - int noConstr=0; - int zeroAtt=0, zeroInstr=0, zeroGlob=0, wrFilebin=0; + long lastStar = -1; + long firstStar = -1; + long lastStarConstr = -1; + long firstStarConstr = -1; + int seqStar = 1; + int withFile = 1; + int noConstr = 0; + int zeroAtt = 0, zeroInstr = 0, zeroGlob = 0, wrFilebin = 0; int constraintFound[MAX_CONSTR][2]; // first: number of file of the constraint; second: row in file - int rowInFile=0; - int noIter=0; - int noCPR=0; - int extConstraint=0, nEqExtConstr=0; - int barConstraint=0, nEqBarConstr=0; - int startingAttColExtConstr,endingAttColExtConstr,numOfExtAttCol,starOverlap=0,numOfExtStar,numOfBarStar; + int rowInFile = 0; + int noIter = 0; + int noCPR = 0; + int extConstraint = 0, nEqExtConstr = 0; + int barConstraint = 0, nEqBarConstr = 0; + int startingAttColExtConstr, endingAttColExtConstr, numOfExtAttCol, starOverlap = 0, numOfExtStar, numOfBarStar; char constranitFoundFileName[MAX_CONSTR][256]; - struct dirent **namelistMI; + struct dirent **namelistMI; struct nullSpace nullSpaceCk; double nullSpaceAttfact; - int autoRun=0; - float memGlobal=0; + int autoRun = 0; + float memGlobal = 0; float memGB; + int nthreads=1; //**QUI modificato a 1 da 0 + int ntasks=0; + /////////////// - for (int i=1;i 0) + timeCPR = testTime; } - if(strcmp (argv[i],"-timeCPR") == 0) + if (strcmp(argv[i], "-timelimit") == 0) { - testTime=atoi(argv[i+1]); - if(testTime >0) timeCPR=testTime; + testTime = atoi(argv[i + 1]); + if (testTime > 0) + timeLimit = testTime; } - if(strcmp (argv[i],"-timelimit") == 0) + + if (strcmp(argv[i], "-itnCPR") == 0) { - testTime=atoi(argv[i+1]); - if(testTime >0) timeLimit=testTime; + testTime = atoi(argv[i + 1]); + if (testTime > 0) + itnCPR = testTime; } - - if(strcmp (argv[i],"-itnCPR") == 0) + if (strcmp(argv[i], "-noCPR") == 0) { - testTime=atoi(argv[i+1]); - if(testTime >0) itnCPR=testTime; + noCPR = 1; } - if(strcmp (argv[i],"-noCPR") == 0) + if (strcmp(argv[i], "-itnlimit") == 0) { - noCPR=1; + testTime = atoi(argv[i + 1]); + if (testTime > 0) + itnLimit = testTime; } - if(strcmp (argv[i],"-itnlimit") == 0) + if(strcmp (argv[i],"-tasks") == 0) { - testTime=atoi(argv[i+1]); - if(testTime >0) itnLimit=testTime; + ntasks=atoi(argv[i+1]); + if(ntasks <=0) ntasks=0; } - if(strcmp (argv[i],"-atol") == 0) + if (strcmp(argv[i], "-atol") == 0) { - double dummy=atof(argv[i+1]); - if(dummy >=0) aTol=dummy; + double dummy = atof(argv[i + 1]); + if (dummy >= 0) + aTol = dummy; } - if(strcmp (argv[i],"-zeroAtt") == 0) + if (strcmp(argv[i], "-zeroAtt") == 0) { - zeroAtt=1; + zeroAtt = 1; } - if(strcmp (argv[i],"-zeroInstr") == 0) + if (strcmp(argv[i], "-zeroInstr") == 0) { - zeroInstr=1; + zeroInstr = 1; } - if(strcmp (argv[i],"-zeroGlob") == 0) + if (strcmp(argv[i], "-zeroGlob") == 0) { - zeroGlob=1; + zeroGlob = 1; } - if(strcmp (argv[i],"-wgic") == 0) + if (strcmp(argv[i], "-wgic") == 0) { - wgInstrCoeff=atoi(argv[i+1]); - if(wgInstrCoeff<=0) wgInstrCoeff=1; + wgInstrCoeff = atoi(argv[i + 1]); + if (wgInstrCoeff <= 0) + wgInstrCoeff = 1; } // inputDir e outputDir - - if(strcmp (argv[i],"-inputDir") == 0) + + if (strcmp(argv[i], "-inputDir") == 0) { - sprintf(inputDir, "%s", argv[i+1]); - inputDirOpt=1; + sprintf(inputDir, "%s", argv[i + 1]); + inputDirOpt = 1; } - - if(strcmp (argv[i],"-outputDir") == 0) + + if (strcmp(argv[i], "-outputDir") == 0) { - sprintf(outputDir, "%s", argv[i+1]); - outputDirOpt=1; + sprintf(outputDir, "%s", argv[i + 1]); + outputDirOpt = 1; } - + // Fine input e output Dir - - if(strcmp (argv[i],"-debug") == 0) - debugMode=1; - - if(strcmp (argv[i],"-wrsol") == 0){ - wrsol=1; - } - - if(strcmp (argv[i],"-noiter") == 0){ - noIter=1; - } - - if(strcmp (argv[i],"-wrfilebin") == 0) - { - sprintf(wrfileDir, "%s", argv[i+1]); - wrFilebin=1; - } - if(strcmp (argv[i],"-extConstr") == 0) - { - extConstraint=1; //extendet external constraint - nEqExtConstr=DEFAULT_EXTCONSTROWS; - extConstrW=atof(argv[i+1]); - if(extConstrW==0.0){ - printf("ERROR: PE=%d -extConstr option given with no value or zero value ==> %le. Aborting\n",myid,extConstrW); - MPI_Abort(MPI_COMM_WORLD,1); + + if (strcmp(argv[i], "-debug") == 0) + debugMode = 1; + + if (strcmp(argv[i], "-wrsol") == 0) + { + wrsol = 1; + } + + if (strcmp(argv[i], "-noiter") == 0) + { + noIter = 1; + } + + if (strcmp(argv[i], "-wrfilebin") == 0) + { + sprintf(wrfileDir, "%s", argv[i + 1]); + wrFilebin = 1; + } + if (strcmp(argv[i], "-extConstr") == 0) + { + extConstraint = 1; //extendet external constraint + nEqExtConstr = DEFAULT_EXTCONSTROWS; + extConstrW = atof(argv[i + 1]); + if (extConstrW == 0.0) + { + printf("ERROR: PE=%d -extConstr option given with no value or zero value ==> %le. Aborting\n", myid, extConstrW); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } } - if(strcmp (argv[i],"-barConstr") == 0) + if (strcmp(argv[i], "-barConstr") == 0) { - barConstraint=1; //extendet external constraint - nEqBarConstr=DEFAULT_BARCONSTROWS; - barConstrW=atof(argv[i+1]); - if(barConstrW==0.0){ - printf("ERROR: PE=%d -barConstr option given with no value or zero value ==> %le. Aborting\n",myid,barConstrW); - MPI_Abort(MPI_COMM_WORLD,1); + barConstraint = 1; //extendet external constraint + nEqBarConstr = DEFAULT_BARCONSTROWS; + barConstrW = atof(argv[i + 1]); + if (barConstrW == 0.0) + { + printf("ERROR: PE=%d -barConstr option given with no value or zero value ==> %le. Aborting\n", myid, barConstrW); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } } - - } - if (extConstraint) noConstr=2; - if (barConstraint) noConstr=2; - - - if(myid==0) - { - printf("Execution of solvergaia Simulator version %s on %d mpi-tasks\n solvergaiaSim ",VER,nproc); - if(idtest)printf ("-IDtest %le ", srIDtest); - if(precond)printf ("-Precond on"); - else printf ("-Precond off"); - if(nfileProc!=3) printf ("-numFilexproc %d ",nfileProc); - if(wrFilebin) printf(" -wrfilebin %s ",wrfileDir); - if(!withFile) printf(" -noFile "); - if(itnLimit>0) printf(" -itnlimit %d ",itnLimit); - if(noCPR>0) printf(" -noCPR "); - if(itnCPR!= DEFAULT_ITNCPR) printf(" -itnCPR %d ", itnCPR); - if(zeroAtt) printf(" -zeroAtt "); - if(noConstr==1) printf(" -noConstr "); - if(zeroInstr) printf(" -zeroInstr "); - if(zeroGlob) printf(" -zeroGlob "); - if(inputDirOpt) printf(" -inputDir %s ", inputDir); - if(outputDirOpt) printf(" -outputDir %s ", outputDir); - if(wrsol) printf(" -wrsol "); - if(noIter) printf(" -noiter "); - if(debugMode) printf(" -debug "); - if(extConstraint)printf("-extConstr %le ",extConstrW); - if(barConstraint)printf("-barConstr %le ",barConstrW); - printf("-wgic %d", wgInstrCoeff); - if(!autoRun) printf(" %s\n",argv[argc-1]); - if(extConstraint && barConstraint) { + } + if (extConstraint) + noConstr = 2; + if (barConstraint) + noConstr = 2; + + if (myid == 0) + { + printf("Execution of solvergaia Simulator version %s on %d mpi-tasks\n solvergaiaSim ", VER, nproc); + if (idtest) + printf("-IDtest %le ", srIDtest); + if (precond) + printf("-Precond on"); + else + printf("-Precond off"); + if (nfileProc != 3) + printf("-numFilexproc %d ", nfileProc); + if (wrFilebin) + printf(" -wrfilebin %s ", wrfileDir); + if (!withFile) + printf(" -noFile "); + if (itnLimit > 0) + printf(" -itnlimit %d ", itnLimit); + if (noCPR > 0) + printf(" -noCPR "); + if (itnCPR != DEFAULT_ITNCPR) + printf(" -itnCPR %d ", itnCPR); + if (zeroAtt) + printf(" -zeroAtt "); + if (noConstr == 1) + printf(" -noConstr "); + if (zeroInstr) + printf(" -zeroInstr "); + if (zeroGlob) + printf(" -zeroGlob "); + if (inputDirOpt) + printf(" -inputDir %s ", inputDir); + if (outputDirOpt) + printf(" -outputDir %s ", outputDir); + if (wrsol) + printf(" -wrsol "); + if (noIter) + printf(" -noiter "); + if (debugMode) + printf(" -debug "); + if(ntasks>0) + printf(" -tasks %d ",ntasks); + if (extConstraint) + printf("-extConstr %le ", extConstrW); + if (barConstraint) + printf("-barConstr %le ", barConstrW); + printf("-wgic %d", wgInstrCoeff); + if (!autoRun) + printf(" %s\n", argv[argc - 1]); + if (extConstraint && barConstraint) + { printf("Error: baricentric anld null space constraints are mutually exclusive. Aborting\n"); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(wrFilebin && strcmp(inputDir, wrfileDir)==0){ + if (wrFilebin && strcmp(inputDir, wrfileDir) == 0) + { printf("inputDir and wrfilebinDir are the same. Execution Aborted.\n"); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); - } + printf ("\n"); } ///////////////// - getcwd(wpath,sizePath); - if(!inputDirOpt) + getcwd(wpath, sizePath); + if (!inputDirOpt) strcpy(inputDir, wpath); - printf("Process %d running on %s\n",myid,processor_name); -#ifdef OMP -#pragma omp parallel + printf("Process %d running on %s\n", myid, processor_name); + ////#ifdef OMP + ///#pragma omp task { - int tid = omp_get_thread_num(); - int nthreads = omp_get_num_threads(); - printf("PE=%d on processor %s total num of threads =%d ID thread =%d \n",myid, processor_name,nthreads,tid); - - } + // int tid = omp_get_thread_num(); + /* + const char* s = getenv("NX_SMP_WORKERS");//? "":"1"; + if(s ==NULL) + s="1"; + nthreads = atoi(s); //QUI** verifica che la chiamata ritorni il numero di tasks allocati!! + */ + + //nthreads sono attuatori + // ntasks sono la mia suddivisione +#ifdef OMP + nthreads=omp_get_num_threads(); //**QUI immagino che nthreads>=1 (se 0 non fiunziona!) #endif - nullSpaceAttfact=1.0*attExtConstrFact*extConstrW; - comlsqr.extConstrW=extConstrW; - comlsqr.nullSpaceAttfact=nullSpaceAttfact; - comlsqr.barConstrW=barConstrW; - - - if(inputDirOpt) - { - if(myid==0) printf("Checking input directory %s ... ", inputDir); - if(!(chdir(inputDir)==0)) { + if(ntasks==0) + ntasks=nthreads; + printf("PE=%d on processor %s total num of threads =%d, ntasks=%d\n",myid, processor_name,nthreads,ntasks); + } + ////#pragma omp taskwait + ////#endif + nullSpaceAttfact = 1.0 * attExtConstrFact * extConstrW; + comlsqr.extConstrW = extConstrW; + comlsqr.nullSpaceAttfact = nullSpaceAttfact; + comlsqr.barConstrW = barConstrW; + + if (inputDirOpt) + { + if (myid == 0) + printf("Checking input directory %s ... ", inputDir); + if (!(chdir(inputDir) == 0)) + { printf("Input directory does not exist. Aborting\n"); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } chdir(wpath); - if(myid==0) printf("success.\n"); + if (myid == 0) + printf("success.\n"); } - inputDirLen=strlen(inputDir); - sprintf(filenameSolProps, "%s", argv[argc-1]); - sprintf(filenameSolPropsFinal,"GsrFinal_%s", argv[argc-1]); - - if(outputDirOpt) + inputDirLen = strlen(inputDir); + sprintf(filenameSolProps, "%s", argv[argc - 1]); + sprintf(filenameSolPropsFinal, "GsrFinal_%s", argv[argc - 1]); + + if (outputDirOpt) { - if(myid==0) printf("Checking output directory %s ...", outputDir); - if(!(chdir(outputDir)==0)) { + if (myid == 0) + printf("Checking output directory %s ...", outputDir); + if (!(chdir(outputDir) == 0)) + { printf("Output directory does not exist on PE %d. Aborting\n", myid); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(myid==0) printf("success.\n"); + if (myid == 0) + printf("success.\n"); } - + ///////////// Initialize the output filename - sprintf(filenameAstroResults, "%s", "GsrAstroParamSolution.bin"); // file storing the Astrometric Parameters - sprintf(filenameAttResults, "%s","GsrAttitudeParamSolution.bin"); // file storing the Attitude Parameters - sprintf(filenameInstrResults, "%s","GsrInstrParamSolution.bin"); // file storing the Instrument Parameters - sprintf(filenameGlobalResults, "%s","GsrGlobalParamSolution.bin"); // file storing the Global Parameters - - -// instrConst=(int *) calloc(DEFAULT_NINSTRINDEXES+1 , sizeof(int)); - instrConst=(int *) calloc(DEFAULT_NINSTRINDEXES , sizeof(int)); - nAttParAxis=4; + sprintf(filenameAstroResults, "%s", "GsrAstroParamSolution.bin"); // file storing the Astrometric Parameters + sprintf(filenameAttResults, "%s", "GsrAttitudeParamSolution.bin"); // file storing the Attitude Parameters + sprintf(filenameInstrResults, "%s", "GsrInstrParamSolution.bin"); // file storing the Instrument Parameters + sprintf(filenameGlobalResults, "%s", "GsrGlobalParamSolution.bin"); // file storing the Global Parameters + + // instrConst=(int *) calloc(DEFAULT_NINSTRINDEXES+1 , sizeof(int)); + instrConst = (int *)calloc(DEFAULT_NINSTRINDEXES, sizeof(int)); + nAttParAxis = 4; //START READING filenameSolProps //START READING filenameSolProps GsrSolProps.dat - if(myid==0 && wrFilebin){ + if (myid == 0 && wrFilebin) + { chdir(wpath); chdir(wrfileDir); - fpSolPropsFileBin=fopen(filenameSolProps,"w"); + fpSolPropsFileBin = fopen(filenameSolProps, "w"); } chdir(wpath); - if(inputDirOpt) chdir(inputDir); - if(myid==0) - { - if (autoRun){ - sphereId=0; - atol= 0.000000; - btol= 0.000000; - conlim= 10000000000000.000000; - itnlim= 2000; - if(itnLimit>0) - itnlim=itnLimit; - damp= 0.000000; - nStar= 10000; - nAstroP= 5; - nAstroPSolved= 5; - nConstrLat= 0; - nConstrLong= 0; - nConstrMuLat= 0; - nConstrMuLong= 0; - nDegFreedomAtt= 230489; - nAttAxes= 3; - instrConst[0]= 1; - instrConst[1]= 62; - instrConst[2]= 1966; - instrConst[3]= 22; - nInstrP= 6; - nInstrPSolved= 6; - lsInstrFlag= 1; - ssInstrFlag= 1; - nuInstrFlag= 1; - maInstrFlag= 1; - nOfInstrConstr= -1; - nElemIC= -1; - nGlobP= 0; - nobs= 2400000; + if (inputDirOpt) + chdir(inputDir); + if (myid == 0) + { + if (autoRun) + { + sphereId = 0; + atol = 0.000000; + btol = 0.000000; + conlim = 10000000000000.000000; + itnlim = 2000; + if (itnLimit > 0) + itnlim = itnLimit; + damp = 0.000000; + nStar = 10000; + nAstroP = 5; + nAstroPSolved = 5; + nConstrLat = 0; + nConstrLong = 0; + nConstrMuLat = 0; + nConstrMuLong = 0; + nDegFreedomAtt = 230489; + nAttAxes = 3; + instrConst[0] = 1; + instrConst[1] = 62; + instrConst[2] = 1966; + instrConst[3] = 22; + nInstrP = 6; + nInstrPSolved = 6; + lsInstrFlag = 1; + ssInstrFlag = 1; + nuInstrFlag = 1; + maInstrFlag = 1; + nOfInstrConstr = -1; + //edited FV + nElemIC = -1; + //nElemIC = 0; + nGlobP = 0; + nobs = 2400000; nAstroParam = nStar * nAstroPSolved; nAttParam = nDegFreedomAtt * nAttAxes; - if(nDegFreedomAtt<4) nAttParAxis=nDegFreedomAtt; - if(nDegFreedomAtt) nAttP = nAttAxes * nAttParAxis; - long nFoVs=1+instrConst[0]; - long nCCDs=instrConst[1]; - long nPixelColumns=instrConst[2]; - long nTimeIntervals=instrConst[3]; + if (nDegFreedomAtt < 4) + nAttParAxis = nDegFreedomAtt; + if (nDegFreedomAtt) + nAttP = nAttAxes * nAttParAxis; + long nFoVs = 1 + instrConst[0]; + long nCCDs = instrConst[1]; + long nPixelColumns = instrConst[2]; + long nTimeIntervals = instrConst[3]; // tot. instr. param. = nCCDs (Cmag) + nFoVs*nCCDs (Cnu) + nCCDs*nPixelColumns (delta_eta) + 3*nFoVs*nCCDs*nTimeIntervals (Delta_eta) + nCCDs*nPixelColumns (delta_zeta) + 3*nFoVs*nCCDs*nTimeIntervals (Delta_zeta) // added flags switch on and off the appropriate kind of parameters - nInstrParam = maInstrFlag*nCCDs+nuInstrFlag*nFoVs*nCCDs+ssInstrFlag*2*nCCDs*nPixelColumns+lsInstrFlag*2*3*nFoVs*nCCDs*nTimeIntervals; - nInstrParamTot = nCCDs+ nFoVs*nCCDs+ 2*nCCDs*nPixelColumns+2*3*nFoVs*nCCDs*nTimeIntervals; - if(nInstrP==0) nInstrParamTot=0; + nInstrParam = maInstrFlag * nCCDs + nuInstrFlag * nFoVs * nCCDs + ssInstrFlag * 2 * nCCDs * nPixelColumns + lsInstrFlag * 2 * 3 * nFoVs * nCCDs * nTimeIntervals; + nInstrParamTot = nCCDs + nFoVs * nCCDs + 2 * nCCDs * nPixelColumns + 2 * 3 * nFoVs * nCCDs * nTimeIntervals; + if (nInstrP == 0) + nInstrParamTot = 0; nGlobalParam = nGlobP; - + nparam = nAstroPSolved + nAttP + nInstrPSolved + nGlobP; - if(memGlobal!=0){ - memGB=simfullram(nStar,nobs,memGlobal,nparam,nAttParam,nInstrParam); - printf("Running with memory %f GB, nStar=%d nobs=%ld\n",memGB, nStar,nobs); + if (memGlobal != 0) + { + memGB = simfullram(nStar, nobs, memGlobal, nparam, nAttParam, nInstrParam); + printf("Running with memory %f GB, nStar=%d nobs=%ld\n", memGB, nStar, nobs); } printf("sphereId= %7ld\n", sphereId); printf("atol= %18.15lf\n", atol); @@ -557,7 +630,7 @@ int main(int argc, char **argv) { printf("nConstrMuLong= %5ld\n", nConstrMuLong); printf("nDegFreedomAtt= %7ld\n", nDegFreedomAtt); printf("nAttAxes= %7hd\n", nAttAxes); - printf("nFovs= %7d ", instrConst[0]+1); + printf("nFovs= %7d ", instrConst[0] + 1); printf("(instrConst[0])= %7d\n", instrConst[0]); printf("nCCDs= %7d\n", instrConst[1]); printf("nPixelColumns= %7d\n", instrConst[2]); @@ -572,616 +645,751 @@ int main(int argc, char **argv) { printf("nElemIC= %7d\n", nElemIC); printf("nGlobP= %7hd\n", nGlobP); printf("nObs= %10ld\n", nobs); - if(wrFilebin){ - fprintf(fpSolPropsFileBin,"sphereId= %ld\n", sphereId); - fprintf(fpSolPropsFileBin,"atol= %lf\n", atol); - fprintf(fpSolPropsFileBin,"btol= %lf\n", btol); - fprintf(fpSolPropsFileBin,"conlim= %lf\n", conlim); - fprintf(fpSolPropsFileBin,"itnlim= %ld\n", itnlim); - fprintf(fpSolPropsFileBin,"damp= %lf\n", damp); - fprintf(fpSolPropsFileBin,"nStar= %ld\n", nStar); - fprintf(fpSolPropsFileBin,"nAstroP= %hd\n", nAstroP); - fprintf(fpSolPropsFileBin,"nAstroPSolved= %hd\n", nAstroPSolved); - fprintf(fpSolPropsFileBin,"nConstrLat= 0\n"); - fprintf(fpSolPropsFileBin,"nConstrLong= 0\n"); - fprintf(fpSolPropsFileBin,"nConstrMuLat= 0\n"); - fprintf(fpSolPropsFileBin,"nConstrMuLong= 0\n"); - fprintf(fpSolPropsFileBin,"nDegFreedomAtt= %ld\n", nDegFreedomAtt); - fprintf(fpSolPropsFileBin,"nAttAxes= %hd\n", nAttAxes); - fprintf(fpSolPropsFileBin,"nFoVs= %d\n", instrConst[0]); - fprintf(fpSolPropsFileBin,"nCCDs= %d\n", instrConst[1]); - fprintf(fpSolPropsFileBin,"nPixelColumns= %d\n", instrConst[2]); - fprintf(fpSolPropsFileBin,"nTimeIntervals= %d\n", instrConst[3]); - fprintf(fpSolPropsFileBin,"nInstrP= %hd\n", nInstrP); - fprintf(fpSolPropsFileBin,"nInstrPSolved= %hd\n", nInstrPSolved); - fprintf(fpSolPropsFileBin,"lsInstrFlag= %hd\n", lsInstrFlag); - fprintf(fpSolPropsFileBin,"ssInstrFlag= %hd\n", ssInstrFlag); - fprintf(fpSolPropsFileBin,"nuInstrFlag= %hd\n", nuInstrFlag); - fprintf(fpSolPropsFileBin,"maInstrFlag= %hd\n", maInstrFlag); - fprintf(fpSolPropsFileBin,"nOfInstrConstr= %d\n", nOfInstrConstr); - fprintf(fpSolPropsFileBin,"nElemIC= %d\n", nElemIC); - fprintf(fpSolPropsFileBin,"nGlobP= %hd\n", nGlobP); - fprintf(fpSolPropsFileBin,"nObs= %ld\n", nobs); + if (wrFilebin) + { + fprintf(fpSolPropsFileBin, "sphereId= %ld\n", sphereId); + fprintf(fpSolPropsFileBin, "atol= %lf\n", atol); + fprintf(fpSolPropsFileBin, "btol= %lf\n", btol); + fprintf(fpSolPropsFileBin, "conlim= %lf\n", conlim); + fprintf(fpSolPropsFileBin, "itnlim= %ld\n", itnlim); + fprintf(fpSolPropsFileBin, "damp= %lf\n", damp); + fprintf(fpSolPropsFileBin, "nStar= %ld\n", nStar); + fprintf(fpSolPropsFileBin, "nAstroP= %hd\n", nAstroP); + fprintf(fpSolPropsFileBin, "nAstroPSolved= %hd\n", nAstroPSolved); + fprintf(fpSolPropsFileBin, "nConstrLat= 0\n"); + fprintf(fpSolPropsFileBin, "nConstrLong= 0\n"); + fprintf(fpSolPropsFileBin, "nConstrMuLat= 0\n"); + fprintf(fpSolPropsFileBin, "nConstrMuLong= 0\n"); + fprintf(fpSolPropsFileBin, "nDegFreedomAtt= %ld\n", nDegFreedomAtt); + fprintf(fpSolPropsFileBin, "nAttAxes= %hd\n", nAttAxes); + fprintf(fpSolPropsFileBin, "nFoVs= %d\n", instrConst[0]); + fprintf(fpSolPropsFileBin, "nCCDs= %d\n", instrConst[1]); + fprintf(fpSolPropsFileBin, "nPixelColumns= %d\n", instrConst[2]); + fprintf(fpSolPropsFileBin, "nTimeIntervals= %d\n", instrConst[3]); + fprintf(fpSolPropsFileBin, "nInstrP= %hd\n", nInstrP); + fprintf(fpSolPropsFileBin, "nInstrPSolved= %hd\n", nInstrPSolved); + fprintf(fpSolPropsFileBin, "lsInstrFlag= %hd\n", lsInstrFlag); + fprintf(fpSolPropsFileBin, "ssInstrFlag= %hd\n", ssInstrFlag); + fprintf(fpSolPropsFileBin, "nuInstrFlag= %hd\n", nuInstrFlag); + fprintf(fpSolPropsFileBin, "maInstrFlag= %hd\n", maInstrFlag); + fprintf(fpSolPropsFileBin, "nOfInstrConstr= %d\n", nOfInstrConstr); + fprintf(fpSolPropsFileBin, "nElemIC= %d\n", nElemIC); + fprintf(fpSolPropsFileBin, "nGlobP= %hd\n", nGlobP); + fprintf(fpSolPropsFileBin, "nObs= %ld\n", nobs); fclose(fpSolPropsFileBin); } - }else{ - fpSolProps=fopen(filenameSolProps,"r"); - if (!fpSolProps) - { - printf("Error while open %s\n",filenameSolProps); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - - if(fscanf(fpSolProps, "sphereId= %ld\n",&sphereId) != 1) { - printf("Error reading sphereId %ld \n",sphereId); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("sphereId= %7ld\n", sphereId); - if(wrFilebin) fprintf(fpSolPropsFileBin,"sphereId= %ld\n", sphereId); - - if(fscanf(fpSolProps, "atol= %lf\n",&atol) != 1) { - printf("Error reading atol %lf \n",atol); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(aTol >=0) - atol=aTol; - printf("atol= %18.15lf\n", atol); - if(wrFilebin) fprintf(fpSolPropsFileBin,"atol= %lf\n", atol); - - if(fscanf(fpSolProps, "btol= %lf\n",&btol) != 1) { - printf("Error reading btol %lf \n",btol); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("btol= %18.15lf\n", btol); - if(wrFilebin) fprintf(fpSolPropsFileBin,"btol= %lf\n", btol); - - if(fscanf(fpSolProps, "conlim= %lf\n",&conlim) != 1) { - printf("Error reading conlim %lf \n",conlim); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("conlim= %18.15le\n", conlim); - if(wrFilebin) fprintf(fpSolPropsFileBin,"conlim= %lf\n", conlim); - - if(fscanf(fpSolProps, "itnlim= %ld\n",&itnlim) != 1) { - printf("Error reading itnlim %ld \n",itnlim); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(itnLimit>0) - itnlim=itnLimit; - printf("itnlim= %7ld\n", itnlim); - if(wrFilebin) fprintf(fpSolPropsFileBin,"itnlim= %ld\n", itnlim); - - if(fscanf(fpSolProps, "damp= %lf\n",&damp) != 1) { - printf("Error reading damp %lf \n",damp); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("damp= %18.15lf\n", damp); - if(wrFilebin) fprintf(fpSolPropsFileBin,"damp= %lf\n", damp); - - if(fscanf(fpSolProps, "nStar= %ld\n",&nStar) != 1) { - printf("Error reading nStar %ld \n",nStar); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nStar= %7ld\n", nStar); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nStar= %ld\n", nStar); - - if(fscanf(fpSolProps, "nAstroP= %hd\n",&nAstroP) != 1) { - printf("Error reading nAstroP %hd \n",nAstroP); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nAstroP= %7hd\n", nAstroP); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nAstroP= %hd\n", nAstroP); - - if(fscanf(fpSolProps, "nAstroPSolved= %hd\n",&nAstroPSolved) != 1) { - printf("Error reading nAstroPSolved %hd \n",nAstroPSolved); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nAstroPSolved= %hd\n", nAstroPSolved); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nAstroPSolved= %hd\n", nAstroPSolved); - - if(fscanf(fpSolProps, "nConstrLat= %ld\n",&nConstrLat) != 1) { - printf("Error reading nConstrLat %ld \n",nConstrLat); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nConstrLat= %5ld\n", nConstrLat); - if(wrFilebin){ - if(noConstr) fprintf(fpSolPropsFileBin,"nConstrLat= 0\n"); - else fprintf(fpSolPropsFileBin,"nConstrLat= %ld\n", nConstrLat); } - if(nConstrLat>0) + else { - constrLatId = (long *) calloc(nConstrLat, sizeof(long)); - for(i=0;i= 0) + atol = aTol; + printf("atol= %18.15lf\n", atol); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "atol= %lf\n", atol); + + if (fscanf(fpSolProps, "btol= %lf\n", &btol) != 1) + { + printf("Error reading btol %lf \n", btol); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(debugMode) printf("constrLatW[nConstrLat-1]=%18.15lf\n",constrLatW[nConstrLat-1]); - if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf\n",constrLatW[nConstrLat-1]); - } - - if(fscanf(fpSolProps, "nConstrLong= %ld\n",&nConstrLong) != 1) { - printf("Error reading nConstrLong %ld \n",nConstrLong); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nConstrLong= %5ld\n", nConstrLong); - if(wrFilebin){ - if(noConstr) fprintf(fpSolPropsFileBin,"nConstrLong= 0\n"); - else fprintf(fpSolPropsFileBin,"nConstrLong= %ld\n", nConstrLong); - } - if(nConstrLong>0) - { - constrLongId = (long *) calloc(nConstrLong, sizeof(long)); - for(i=0;i 0) + itnlim = itnLimit; + printf("itnlim= %7ld\n", itnlim); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "itnlim= %ld\n", itnlim); + + if (fscanf(fpSolProps, "damp= %lf\n", &damp) != 1) + { + printf("Error reading damp %lf \n", damp); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); } - if(fscanf(fpSolProps, "%lf\n",&constrLongW[nConstrLong-1]) != 1) { - printf("Error reading constrLongW[nConstrLong-1] %lf \n", constrLongW[nConstrLong-1]); - MPI_Abort(MPI_COMM_WORLD,1); + printf("damp= %18.15lf\n", damp); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "damp= %lf\n", damp); + + if (fscanf(fpSolProps, "nStar= %ld\n", &nStar) != 1) + { + printf("Error reading nStar %ld \n", nStar); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(debugMode) printf("constrLongW[nConstrLong-1]=%18.15lf\n",constrLongW[nConstrLong-1]); - if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf\n",constrLongW[nConstrLong-1]); - } - - if(fscanf(fpSolProps, "nConstrMuLat= %ld\n",&nConstrMuLat) != 1) { - printf("Error reading nConstrMuLat %ld \n",nConstrMuLat); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nConstrMuLat= %5ld\n", nConstrMuLat); - if(wrFilebin){ - if(noConstr) fprintf(fpSolPropsFileBin,"nConstrMuLat= 0\n"); - else fprintf(fpSolPropsFileBin,"nConstrMuLat= %ld\n", nConstrMuLat); - } - - if(nConstrMuLat>0) - { - constrMuLatId = (long *) calloc(nConstrMuLat, sizeof(long)); - for(i=0;i 0) + { + constrLatId = (long *)calloc(nConstrLat, sizeof(long)); + for (i = 0; i < nConstrLat - 1; i++) + { + if (fscanf(fpSolProps, "%ld ", &constrLatId[i]) != 1) + { + printf("Error reading constrLatId[%d] %ld \n", i, constrLatId[i]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrLatId[%d]=%7ld ", i, constrLatId[i]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%ld ", constrLatId[i]); + } + if (fscanf(fpSolProps, "%ld\n", &constrLatId[nConstrLat - 1]) != 1) + { + printf("Error reading constrLatId[nConstrLat-1] %ld \n", constrLatId[nConstrLat - 1]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrLatId[nConstrLat-1]=%7ld\n", constrLatId[nConstrLat - 1]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%ld\n", constrLatId[nConstrLat - 1]); + + constrLatW = (double *)calloc(nConstrLat, sizeof(double)); + for (i = 0; i < nConstrLat - 1; i++) + { + if (fscanf(fpSolProps, "%lf ", &constrLatW[i]) != 1) + { + printf("Error reading constrLatW[%d] %lf \n", i, constrLatW[i]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrLatW[%d]=%18.15le ", i, constrLatW[i]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%lf ", constrLatW[i]); + } + if (fscanf(fpSolProps, "%lf\n", &constrLatW[nConstrLat - 1]) != 1) + { + printf("Error reading constrLatW[nConstrLat-1] %lf \n", constrLatW[nConstrLat - 1]); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(debugMode) printf("constrMuLatW[%d]=%18.15le ",i,constrMuLatW[i]); - if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf ",constrMuLatW[i]); + if (debugMode) + printf("constrLatW[nConstrLat-1]=%18.15lf\n", constrLatW[nConstrLat - 1]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%lf\n", constrLatW[nConstrLat - 1]); } - if(fscanf(fpSolProps, "%lf\n",&constrMuLatW[nConstrMuLat-1]) != 1) { - printf("Error reading constrMuLatW[nConstrMuLat-1] %lf \n", constrMuLatW[nConstrMuLat-1]); - MPI_Abort(MPI_COMM_WORLD,1); + + if (fscanf(fpSolProps, "nConstrLong= %ld\n", &nConstrLong) != 1) + { + printf("Error reading nConstrLong %ld \n", nConstrLong); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(debugMode) printf("constrMuLatW[nConstrMuLat-1]=%18.15lf\n",constrMuLatW[nConstrMuLat-1]); - if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf\n",constrMuLatW[nConstrMuLat-1]); - - } - - if(fscanf(fpSolProps, "nConstrMuLong= %ld\n",&nConstrMuLong) != 1) { - printf("Error reading nConstrMuLong %ld \n",nConstrMuLong); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nConstrMuLong= %5ld\n", nConstrMuLong); - if(wrFilebin){ - if(noConstr) fprintf(fpSolPropsFileBin,"nConstrMuLong= 0\n"); - else fprintf(fpSolPropsFileBin,"nConstrMuLong= %ld\n", nConstrMuLong); - } - if(nConstrMuLong>0) - { - constrMuLongId = (long *) calloc(nConstrMuLong, sizeof(long)); - for(i=0;i 0) + { + constrLongId = (long *)calloc(nConstrLong, sizeof(long)); + for (i = 0; i < nConstrLong - 1; i++) + { + if (fscanf(fpSolProps, "%ld ", &constrLongId[i]) != 1) + { + printf("Error reading constrLongId[%d] %ld \n", i, constrLongId[i]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrLongId[%d]=%7ld ", i, constrLongId[i]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%ld ", constrLongId[i]); + } + if (fscanf(fpSolProps, "%ld\n", &constrLongId[nConstrLong - 1]) != 1) + { + printf("Error reading constrLongId[nConstrLong-1] %ld \n", constrLongId[nConstrLong - 1]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrLongId[nConstrLong-1]=%7ld\n", constrLongId[nConstrLong - 1]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%ld\n", constrLongId[nConstrLong - 1]); + + constrLongW = (double *)calloc(nConstrLong, sizeof(double)); + for (i = 0; i < nConstrLong - 1; i++) + { + if (fscanf(fpSolProps, "%lf ", &constrLongW[i]) != 1) + { + printf("Error reading constrLongW[%d] %lf \n", i, constrLongW[i]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrLongW[%d]=%18.15le ", i, constrLongW[i]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%lf ", constrLongW[i]); + } + if (fscanf(fpSolProps, "%lf\n", &constrLongW[nConstrLong - 1]) != 1) + { + printf("Error reading constrLongW[nConstrLong-1] %lf \n", constrLongW[nConstrLong - 1]); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(debugMode) printf("constrMuLongId[%d]=%7ld ",i,constrMuLongId[i]); - if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%7ld ",constrMuLongId[i]); + if (debugMode) + printf("constrLongW[nConstrLong-1]=%18.15lf\n", constrLongW[nConstrLong - 1]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%lf\n", constrLongW[nConstrLong - 1]); } - if(fscanf(fpSolProps, "%ld\n",&constrMuLongId[nConstrMuLong-1]) != 1) { - printf("Error reading constrMuLongId[nConstrMuLong-1] %ld \n", constrMuLongId[nConstrMuLong-1]); - MPI_Abort(MPI_COMM_WORLD,1); + + if (fscanf(fpSolProps, "nConstrMuLat= %ld\n", &nConstrMuLat) != 1) + { + printf("Error reading nConstrMuLat %ld \n", nConstrMuLat); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(debugMode) printf("constrMuLongId[nConstrMuLong-1]=%7ld\n",constrMuLongId[nConstrMuLong-1]); - if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%ld\n",constrMuLongId[nConstrMuLong-1]); - - constrMuLongW = (double *) calloc(nConstrMuLong, sizeof(double)); - for(i=0;i 0) + { + constrMuLatId = (long *)calloc(nConstrMuLat, sizeof(long)); + for (i = 0; i < nConstrMuLat - 1; i++) + { + if (fscanf(fpSolProps, "%ld ", &constrMuLatId[i]) != 1) + { + printf("Error reading constrMuLatId[%d] %ld \n", i, constrMuLatId[i]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrMuLatId[%d]=%7ld ", i, constrMuLatId[i]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%ld ", constrMuLatId[i]); + } + if (fscanf(fpSolProps, "%ld\n", &constrMuLatId[nConstrMuLat - 1]) != 1) + { + printf("Error reading constrMuLatId[nConstrMuLat-1] %ld \n", constrMuLatId[nConstrMuLat - 1]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrMuLatId[nConstrMuLat-1]=%7ld\n", constrMuLatId[nConstrMuLat - 1]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%ld\n", constrMuLatId[nConstrMuLat - 1]); + + constrMuLatW = (double *)calloc(nConstrMuLat, sizeof(double)); + for (i = 0; i < nConstrMuLat - 1; i++) + { + if (fscanf(fpSolProps, "%lf ", &constrMuLatW[i]) != 1) + { + printf("Error reading constrMuLatW[%d] %lf \n", i, constrMuLatW[i]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrMuLatW[%d]=%18.15le ", i, constrMuLatW[i]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%lf ", constrMuLatW[i]); + } + if (fscanf(fpSolProps, "%lf\n", &constrMuLatW[nConstrMuLat - 1]) != 1) + { + printf("Error reading constrMuLatW[nConstrMuLat-1] %lf \n", constrMuLatW[nConstrMuLat - 1]); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(debugMode) printf("constrMuLongW[%d]=%18.15le ",i,constrMuLongW[i]); - if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf ",constrMuLongW[i]); + if (debugMode) + printf("constrMuLatW[nConstrMuLat-1]=%18.15lf\n", constrMuLatW[nConstrMuLat - 1]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%lf\n", constrMuLatW[nConstrMuLat - 1]); } - if(fscanf(fpSolProps, "%lf\n",&constrMuLongW[nConstrMuLong-1]) != 1) { - printf("Error reading constrMuLongW[nConstrMuLong-1] %lf \n", constrMuLongW[nConstrMuLong-1]); - MPI_Abort(MPI_COMM_WORLD,1); + + if (fscanf(fpSolProps, "nConstrMuLong= %ld\n", &nConstrMuLong) != 1) + { + printf("Error reading nConstrMuLong %ld \n", nConstrMuLong); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - if(debugMode) printf("constrMuLongW[nConstrMuLong-1]=%18.15lf\n",constrMuLongW[nConstrMuLong-1]); - if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf\n",constrMuLongW[nConstrMuLong-1]); - } - - if(fscanf(fpSolProps, "nDegFreedomAtt= %ld\n",&nDegFreedomAtt) != 1) { - printf("Error reading nDegFreedomAtt %ld \n",nDegFreedomAtt); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(zeroAtt) nDegFreedomAtt=0; - printf("nDegFreedomAtt= %7ld\n", nDegFreedomAtt); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nDegFreedomAtt= %ld\n", nDegFreedomAtt); - - if(fscanf(fpSolProps, "nAttAxes= %hd\n",&nAttAxes) != 1) { - printf("Error reading nAttAxes %hd \n",nAttAxes); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(zeroAtt) nAttAxes=0; - printf("nAttAxes= %7hd\n", nAttAxes); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nAttAxes= %hd\n", nAttAxes); - - if(fscanf(fpSolProps, "nFoVs= %d\n",&instrConst[0]) != 1) { - printf("Error reading nFoVs %d \n",instrConst[0]); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(zeroInstr) instrConst[0]=0; - // instrConst[0] must be 0 or 1 because nFoVs = 2 for Gaia and nFoVs=1+instrConst[0] - if(instrConst[0]<0 || instrConst[0]>1){ - printf("Execution aborted with invalid nFovs=%d\n",instrConst[0]+1); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - - } - printf("nFovs= %7d ", instrConst[0]+1); - printf("(instrConst[0])= %7d\n", instrConst[0]); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nFoVs= %d\n", instrConst[0]); - - if(fscanf(fpSolProps, "nCCDs= %d\n",&instrConst[1]) != 1) { - printf("Error reading nCCDs %d \n",instrConst[1]); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(zeroInstr) instrConst[1]=0; - printf("nCCDs= %7d\n", instrConst[1]); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nCCDs= %d\n", instrConst[1]); - - if(fscanf(fpSolProps, "nPixelColumns= %d\n",&instrConst[2]) != 1) { - printf("Error reading nPixelColumns %d \n",instrConst[2]); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(zeroInstr) instrConst[2]=0; - printf("nPixelColumns= %7d\n", instrConst[2]); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nPixelColumns= %d\n", instrConst[2]); - - if(fscanf(fpSolProps, "nTimeIntervals= %d\n",&instrConst[3]) != 1) { - printf("Error reading nTimeIntervals %d \n",instrConst[3]); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(zeroInstr) instrConst[3]=0; - printf("nTimeIntervals= %7d\n", instrConst[3]); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nTimeIntervals= %d\n", instrConst[3]); - - if(fscanf(fpSolProps, "nInstrP= %hd\n",&nInstrP) != 1) { - printf("Error reading nInstrP %hd \n",nInstrP); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(zeroInstr) nInstrP=0; - printf("nInstrP= %7hd\n", nInstrP); - if(nInstrP!=0 && nInstrP !=6){ - printf("Execution aborted with invalid nInstrP\n"); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(wrFilebin) fprintf(fpSolPropsFileBin,"nInstrP= %hd\n", nInstrP); - if(fscanf(fpSolProps, "nInstrPSolved= %hd\n",&nInstrPSolved) != 1) - { - printf("Error reading nInstrPSolved %hd \n",nInstrPSolved); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nInstrPSolved= %7hd\n", nInstrPSolved); - - if(nInstrPSolved<0 || nInstrPSolved >6){ - printf("Execution aborted with invalid nInstrPSolved\n"); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(wrFilebin) fprintf(fpSolPropsFileBin,"nInstrPSolved= %hd\n", nInstrPSolved); - - if(fscanf(fpSolProps, "lsInstrFlag= %d\n",&lsInstrFlag) != 1) - { - printf("Error reading lsInstrFlag %d \n",lsInstrFlag); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("lsInstrFlag= %7d\n", lsInstrFlag); - - if(lsInstrFlag<0 || lsInstrFlag >1){ - printf("Execution aborted with invalid lsInstrFlag\n"); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(wrFilebin) fprintf(fpSolPropsFileBin,"lsInstrFlag= %hd\n", lsInstrFlag); - if(fscanf(fpSolProps, "ssInstrFlag= %d\n",&ssInstrFlag) != 1) - { - printf("Error reading ssInstrFlag %d \n",ssInstrFlag); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("ssInstrFlag= %7d\n", ssInstrFlag); - - if(ssInstrFlag<0 || ssInstrFlag >1){ - printf("Execution aborted with invalid ssInstrFlag\n"); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(wrFilebin) fprintf(fpSolPropsFileBin,"ssInstrFlag= %hd\n", ssInstrFlag); - if(fscanf(fpSolProps, "nuInstrFlag= %d\n",&nuInstrFlag) != 1) - { - printf("Error reading nuInstrFlag %d \n",nuInstrFlag); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nuInstrFlag= %7d\n", nuInstrFlag); - - if(nuInstrFlag<0 || nuInstrFlag >1){ - printf("Execution aborted with invalid lsInstrFlag\n"); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(wrFilebin) fprintf(fpSolPropsFileBin,"nuInstrFlag= %hd\n", nuInstrFlag); - if(fscanf(fpSolProps, "maInstrFlag= %d\n",&maInstrFlag) != 1) - { - printf("Error reading maInstrFlag %d \n",maInstrFlag); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("maInstrFlag= %7d\n", maInstrFlag); - - if(maInstrFlag<0 || maInstrFlag >1){ - printf("Execution aborted with invalid maInstrFlag\n"); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(wrFilebin) fprintf(fpSolPropsFileBin,"maInstrFlag= %hd\n", maInstrFlag); + printf("nConstrMuLong= %5ld\n", nConstrMuLong); + if (wrFilebin) + { + if (noConstr) + fprintf(fpSolPropsFileBin, "nConstrMuLong= 0\n"); + else + fprintf(fpSolPropsFileBin, "nConstrMuLong= %ld\n", nConstrMuLong); + } + if (nConstrMuLong > 0) + { + constrMuLongId = (long *)calloc(nConstrMuLong, sizeof(long)); + for (i = 0; i < nConstrMuLong - 1; i++) + { + if (fscanf(fpSolProps, "%ld ", &constrMuLongId[i]) != 1) + { + printf("Error reading constrMuLongId[%d] %ld \n", i, constrMuLongId[i]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrMuLongId[%d]=%7ld ", i, constrMuLongId[i]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%7ld ", constrMuLongId[i]); + } + if (fscanf(fpSolProps, "%ld\n", &constrMuLongId[nConstrMuLong - 1]) != 1) + { + printf("Error reading constrMuLongId[nConstrMuLong-1] %ld \n", constrMuLongId[nConstrMuLong - 1]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrMuLongId[nConstrMuLong-1]=%7ld\n", constrMuLongId[nConstrMuLong - 1]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%ld\n", constrMuLongId[nConstrMuLong - 1]); - if(nInstrPSolved != maInstrFlag+nuInstrFlag+ssInstrFlag+3*lsInstrFlag){ - printf("Execution aborted invalid nInstrPSolved=%d. It should be =%d\n",nInstrPSolved,maInstrFlag+nuInstrFlag+ssInstrFlag+3*lsInstrFlag); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - - } - - if(fscanf(fpSolProps, "nOfInstrConstr= %d\n",&nOfInstrConstr) != 1) - { - printf("Error reading nOfInstrConstr %d \n",nOfInstrConstr); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nOfInstrConstr= %7d\n", nOfInstrConstr); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nOfInstrConstr= %d\n", nOfInstrConstr); - - if(fscanf(fpSolProps, "nElemIC= %d\n",&nElemIC) != 1) + constrMuLongW = (double *)calloc(nConstrMuLong, sizeof(double)); + for (i = 0; i < nConstrMuLong - 1; i++) + { + if (fscanf(fpSolProps, "%lf ", &constrMuLongW[i]) != 1) + { + printf("Error reading constrMuLongW[%d] %lf \n", i, constrMuLongW[i]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrMuLongW[%d]=%18.15le ", i, constrMuLongW[i]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%lf ", constrMuLongW[i]); + } + if (fscanf(fpSolProps, "%lf\n", &constrMuLongW[nConstrMuLong - 1]) != 1) + { + printf("Error reading constrMuLongW[nConstrMuLong-1] %lf \n", constrMuLongW[nConstrMuLong - 1]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (debugMode) + printf("constrMuLongW[nConstrMuLong-1]=%18.15lf\n", constrMuLongW[nConstrMuLong - 1]); + if (wrFilebin && !noConstr) + fprintf(fpSolPropsFileBin, "%lf\n", constrMuLongW[nConstrMuLong - 1]); + } + + if (fscanf(fpSolProps, "nDegFreedomAtt= %ld\n", &nDegFreedomAtt) != 1) + { + printf("Error reading nDegFreedomAtt %ld \n", nDegFreedomAtt); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (zeroAtt) + nDegFreedomAtt = 0; + printf("nDegFreedomAtt= %7ld\n", nDegFreedomAtt); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nDegFreedomAtt= %ld\n", nDegFreedomAtt); + + if (fscanf(fpSolProps, "nAttAxes= %hd\n", &nAttAxes) != 1) + { + printf("Error reading nAttAxes %hd \n", nAttAxes); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (zeroAtt) + nAttAxes = 0; + printf("nAttAxes= %7hd\n", nAttAxes); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nAttAxes= %hd\n", nAttAxes); + + if (fscanf(fpSolProps, "nFoVs= %d\n", &instrConst[0]) != 1) + { + printf("Error reading nFoVs %d \n", instrConst[0]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (zeroInstr) + instrConst[0] = 0; + // instrConst[0] must be 0 or 1 because nFoVs = 2 for Gaia and nFoVs=1+instrConst[0] + if (instrConst[0] < 0 || instrConst[0] > 1) + { + printf("Execution aborted with invalid nFovs=%d\n", instrConst[0] + 1); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + printf("nFovs= %7d ", instrConst[0] + 1); + printf("(instrConst[0])= %7d\n", instrConst[0]); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nFoVs= %d\n", instrConst[0]); + + if (fscanf(fpSolProps, "nCCDs= %d\n", &instrConst[1]) != 1) + { + printf("Error reading nCCDs %d \n", instrConst[1]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (zeroInstr) + instrConst[1] = 0; + printf("nCCDs= %7d\n", instrConst[1]); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nCCDs= %d\n", instrConst[1]); + + if (fscanf(fpSolProps, "nPixelColumns= %d\n", &instrConst[2]) != 1) + { + printf("Error reading nPixelColumns %d \n", instrConst[2]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (zeroInstr) + instrConst[2] = 0; + printf("nPixelColumns= %7d\n", instrConst[2]); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nPixelColumns= %d\n", instrConst[2]); + + if (fscanf(fpSolProps, "nTimeIntervals= %d\n", &instrConst[3]) != 1) + { + printf("Error reading nTimeIntervals %d \n", instrConst[3]); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (zeroInstr) + instrConst[3] = 0; + printf("nTimeIntervals= %7d\n", instrConst[3]); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nTimeIntervals= %d\n", instrConst[3]); + + if (fscanf(fpSolProps, "nInstrP= %hd\n", &nInstrP) != 1) + { + printf("Error reading nInstrP %hd \n", nInstrP); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (zeroInstr) + nInstrP = 0; + printf("nInstrP= %7hd\n", nInstrP); + if (nInstrP != 0 && nInstrP != 6) + { + printf("Execution aborted with invalid nInstrP\n"); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nInstrP= %hd\n", nInstrP); + if (fscanf(fpSolProps, "nInstrPSolved= %hd\n", &nInstrPSolved) != 1) + { + printf("Error reading nInstrPSolved %hd \n", nInstrPSolved); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + printf("nInstrPSolved= %7hd\n", nInstrPSolved); + + if (nInstrPSolved < 0 || nInstrPSolved > 6) + { + printf("Execution aborted with invalid nInstrPSolved\n"); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nInstrPSolved= %hd\n", nInstrPSolved); + + if (fscanf(fpSolProps, "lsInstrFlag= %d\n", &lsInstrFlag) != 1) + { + printf("Error reading lsInstrFlag %d \n", lsInstrFlag); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + printf("lsInstrFlag= %7d\n", lsInstrFlag); + + if (lsInstrFlag < 0 || lsInstrFlag > 1) + { + printf("Execution aborted with invalid lsInstrFlag\n"); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (wrFilebin) + fprintf(fpSolPropsFileBin, "lsInstrFlag= %hd\n", lsInstrFlag); + if (fscanf(fpSolProps, "ssInstrFlag= %d\n", &ssInstrFlag) != 1) + { + printf("Error reading ssInstrFlag %d \n", ssInstrFlag); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + printf("ssInstrFlag= %7d\n", ssInstrFlag); + + if (ssInstrFlag < 0 || ssInstrFlag > 1) + { + printf("Execution aborted with invalid ssInstrFlag\n"); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (wrFilebin) + fprintf(fpSolPropsFileBin, "ssInstrFlag= %hd\n", ssInstrFlag); + if (fscanf(fpSolProps, "nuInstrFlag= %d\n", &nuInstrFlag) != 1) + { + printf("Error reading nuInstrFlag %d \n", nuInstrFlag); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + printf("nuInstrFlag= %7d\n", nuInstrFlag); + + if (nuInstrFlag < 0 || nuInstrFlag > 1) + { + printf("Execution aborted with invalid lsInstrFlag\n"); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nuInstrFlag= %hd\n", nuInstrFlag); + if (fscanf(fpSolProps, "maInstrFlag= %d\n", &maInstrFlag) != 1) + { + printf("Error reading maInstrFlag %d \n", maInstrFlag); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + printf("maInstrFlag= %7d\n", maInstrFlag); + + if (maInstrFlag < 0 || maInstrFlag > 1) + { + printf("Execution aborted with invalid maInstrFlag\n"); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (wrFilebin) + fprintf(fpSolPropsFileBin, "maInstrFlag= %hd\n", maInstrFlag); + + if (nInstrPSolved != maInstrFlag + nuInstrFlag + ssInstrFlag + 3 * lsInstrFlag) + { + printf("Execution aborted invalid nInstrPSolved=%d. It should be =%d\n", nInstrPSolved, maInstrFlag + nuInstrFlag + ssInstrFlag + 3 * lsInstrFlag); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + + if (fscanf(fpSolProps, "nOfInstrConstr= %d\n", &nOfInstrConstr) != 1) + { + printf("Error reading nOfInstrConstr %d \n", nOfInstrConstr); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + printf("nOfInstrConstr= %7d\n", nOfInstrConstr); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nOfInstrConstr= %d\n", nOfInstrConstr); + + if (fscanf(fpSolProps, "nElemIC= %d\n", &nElemIC) != 1) + { + printf("Error reading nElemIC %d \n", nElemIC); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + printf("nElemIC= %7d\n", nElemIC); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nElemIC= %d\n", nElemIC); + + if (fscanf(fpSolProps, "nGlobP= %hd\n", &nGlobP) != 1) + { + printf("Error reading nGlobP %hd \n", nGlobP); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + if (zeroGlob) + nGlobP = 0; + printf("nGlobP= %7hd\n", nGlobP); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nGlobP= %hd\n", nGlobP); + + if (fscanf(fpSolProps, "nObs= %ld\n", &nobs) != 1) + { + printf("Error reading nObs %ld \n", nobs); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); + } + printf("nObs= %10ld\n", nobs); + if (wrFilebin) + fprintf(fpSolPropsFileBin, "nObs= %ld\n", nobs); + + fclose(fpSolProps); + if (wrFilebin) + fclose(fpSolPropsFileBin); + } //if (autoRun) else + } + + endTime = MPI_Wtime(); + if ((nDegFreedomAtt == 0 && nAttAxes > 0) || (nDegFreedomAtt > 0 && nAttAxes == 0)) + { + if (myid == 0) { - printf("Error reading nElemIC %d \n",nElemIC); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - printf("nElemIC= %7d\n", nElemIC); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nElemIC= %d\n", nElemIC); - - - if(fscanf(fpSolProps, "nGlobP= %hd\n",&nGlobP) != 1) { - printf("Error reading nGlobP %hd \n",nGlobP); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); - } - if(zeroGlob) nGlobP=0; - printf("nGlobP= %7hd\n", nGlobP); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nGlobP= %hd\n", nGlobP); - - if(fscanf(fpSolProps, "nObs= %ld\n",&nobs) != 1) { - printf("Error reading nObs %ld \n",nobs); - MPI_Abort(MPI_COMM_WORLD,1); + printf("inconsistent values for nDegFreedomAtt=%ld and nAttaxes=%d\n", nDegFreedomAtt, nAttAxes); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - printf("nObs= %10ld\n", nobs); - if(wrFilebin) fprintf(fpSolPropsFileBin,"nObs= %ld\n", nobs); - - fclose(fpSolProps); - if(wrFilebin) fclose(fpSolPropsFileBin); - }//if (autoRun) else - } - - endTime=MPI_Wtime(); - if((nDegFreedomAtt==0 && nAttAxes>0) || (nDegFreedomAtt>0 && nAttAxes==0) ){ - if(myid==0){ - printf("inconsistent values for nDegFreedomAtt=%ld and nAttaxes=%d\n",nDegFreedomAtt,nAttAxes); - MPI_Abort(MPI_COMM_WORLD, 1); - exit(EXIT_FAILURE); - } } - - if(myid==0) printf("Time to read initial parameters =%f sec.\n",endTime-startTime); + + if (myid == 0) + printf("Time to read initial parameters =%f sec.\n", endTime - startTime); // Start section for parameter broadcast - MPI_Bcast( &atol, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast( &btol, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast( &conlim, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast( &damp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast( &nAstroP, 1, MPI_SHORT, 0, MPI_COMM_WORLD); - MPI_Bcast( &nAstroPSolved, 1, MPI_SHORT, 0, MPI_COMM_WORLD); - MPI_Bcast( &nInstrP, 1, MPI_SHORT, 0, MPI_COMM_WORLD); - MPI_Bcast( &nInstrPSolved, 1, MPI_SHORT, 0, MPI_COMM_WORLD); - MPI_Bcast( &lsInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast( &ssInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast( &nuInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast( &maInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast( &nElemIC, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast( &nOfInstrConstr, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast( &nGlobP, 1, MPI_SHORT, 0, MPI_COMM_WORLD); - MPI_Bcast( &itnlim, 1, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( &nStar, 1, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( &nDegFreedomAtt, 1, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( &nAttAxes, 1, MPI_SHORT, 0, MPI_COMM_WORLD); - MPI_Bcast( &instrSetUp, 1, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( instrConst, DEFAULT_NINSTRINDEXES , MPI_INT, 0, MPI_COMM_WORLD); // errore - MPI_Bcast( &nobs, 1, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( &nConstrLat, 1, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( &nConstrLong, 1, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( &nConstrMuLat, 1, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( &nConstrMuLong, 1, MPI_LONG, 0, MPI_COMM_WORLD); - if(myid !=0) - { - constrLatId=(long *) calloc(nConstrLat , sizeof(long)); - constrLatW=(double *) calloc(nConstrLat , sizeof(double)); - constrLongId=(long *) calloc(nConstrLong , sizeof(long)); - constrLongW=(double *) calloc(nConstrLong , sizeof(double)); - constrMuLatId=(long *) calloc(nConstrMuLat , sizeof(long)); - constrMuLatW=(double *) calloc(nConstrMuLat , sizeof(double)); - constrMuLongId=(long *) calloc(nConstrMuLong , sizeof(long)); - constrMuLongW=(double *) calloc(nConstrMuLong , sizeof(double)); - } - MPI_Bcast( constrLatId, nConstrLat, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( constrLatW, nConstrLat, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast( constrLongId, nConstrLong, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( constrLongW, nConstrLong, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast( constrMuLatId, nConstrMuLat, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( constrMuLatW, nConstrMuLat, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast( constrMuLongId, nConstrMuLong, MPI_LONG, 0, MPI_COMM_WORLD); - MPI_Bcast( constrMuLongW, nConstrMuLong, MPI_DOUBLE, 0, MPI_COMM_WORLD); - - if(nInstrPSolved==0) - zeroInstr=1; + MPI_Bcast(&atol, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&btol, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&conlim, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&damp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&nAstroP, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast(&nAstroPSolved, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast(&nInstrP, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast(&nInstrPSolved, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast(&lsInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&ssInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&nuInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&maInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&nElemIC, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&nOfInstrConstr, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&nGlobP, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast(&itnlim, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(&nStar, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(&nDegFreedomAtt, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(&nAttAxes, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast(&instrSetUp, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(instrConst, DEFAULT_NINSTRINDEXES, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&nobs, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(&nConstrLat, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(&nConstrLong, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(&nConstrMuLat, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(&nConstrMuLong, 1, MPI_LONG, 0, MPI_COMM_WORLD); + if (myid != 0) + { + constrLatId = (long *)calloc(nConstrLat, sizeof(long)); + constrLatW = (double *)calloc(nConstrLat, sizeof(double)); + constrLongId = (long *)calloc(nConstrLong, sizeof(long)); + constrLongW = (double *)calloc(nConstrLong, sizeof(double)); + constrMuLatId = (long *)calloc(nConstrMuLat, sizeof(long)); + constrMuLatW = (double *)calloc(nConstrMuLat, sizeof(double)); + constrMuLongId = (long *)calloc(nConstrMuLong, sizeof(long)); + constrMuLongW = (double *)calloc(nConstrMuLong, sizeof(double)); + } + MPI_Bcast(constrLatId, nConstrLat, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(constrLatW, nConstrLat, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(constrLongId, nConstrLong, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(constrLongW, nConstrLong, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(constrMuLatId, nConstrMuLat, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(constrMuLatW, nConstrMuLat, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(constrMuLatW, nConstrMuLat, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(constrMuLongId, nConstrMuLong, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast(constrMuLongW, nConstrMuLong, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + if (nInstrPSolved == 0) + zeroInstr = 1; /////////// Multiplicity of matrixIndex - + int multMI; - - if(nAstroPSolved) - multMI=2; - else{ - multMI=1; - extConstraint=0; - barConstraint=0; + + if (nAstroPSolved) + multMI = 2; + else + { + multMI = 1; + extConstraint = 0; + barConstraint = 0; } nAstroParam = nStar * nAstroPSolved; nAttParam = nDegFreedomAtt * nAttAxes; - if(nDegFreedomAtt<4) nAttParAxis=nDegFreedomAtt; - if(nDegFreedomAtt) nAttP = nAttAxes * nAttParAxis; - long nFoVs=1+instrConst[0]; - long nCCDs=instrConst[1]; - long nPixelColumns=instrConst[2]; - long nTimeIntervals=instrConst[3]; + if (nDegFreedomAtt < 4) + nAttParAxis = nDegFreedomAtt; + if (nDegFreedomAtt) + nAttP = nAttAxes * nAttParAxis; + long nFoVs = 1 + instrConst[0]; + long nCCDs = instrConst[1]; + long nPixelColumns = instrConst[2]; + long nTimeIntervals = instrConst[3]; // tot. instr. param. = nCCDs (Cmag) + nFoVs*nCCDs (Cnu) + nCCDs*nPixelColumns (delta_eta) + 3*nFoVs*nCCDs*nTimeIntervals (Delta_eta) + nCCDs*nPixelColumns (delta_zeta) + 3*nFoVs*nCCDs*nTimeIntervals (Delta_zeta) // added flags switch on and off the appropriate kind of parameters - nInstrParam = maInstrFlag*nCCDs+nuInstrFlag*nFoVs*nCCDs+ssInstrFlag*2*nCCDs*nPixelColumns+lsInstrFlag*2*3*nFoVs*nCCDs*nTimeIntervals; - nInstrParamTot = nCCDs+ nFoVs*nCCDs+ 2*nCCDs*nPixelColumns+2*3*nFoVs*nCCDs*nTimeIntervals; + nInstrParam = maInstrFlag * nCCDs + nuInstrFlag * nFoVs * nCCDs + ssInstrFlag * 2 * nCCDs * nPixelColumns + lsInstrFlag * 2 * 3 * nFoVs * nCCDs * nTimeIntervals; + nInstrParamTot = nCCDs + nFoVs * nCCDs + 2 * nCCDs * nPixelColumns + 2 * 3 * nFoVs * nCCDs * nTimeIntervals; nGlobalParam = nGlobP; - - if(nElemIC<0 || nOfInstrConstr <0){ - - - nOfInstrConstrLSAL = lsInstrFlag*nTimeIntervals; - nElemICLSAL = nFoVs*nCCDs; - nOfInstrConstrLSAC = lsInstrFlag*nFoVs*nTimeIntervals; + + if (nElemIC < 0 || nOfInstrConstr < 0) + { + nOfInstrConstrLSAL = lsInstrFlag * nTimeIntervals; + nElemICLSAL = nFoVs * nCCDs; + nOfInstrConstrLSAC = lsInstrFlag * nFoVs * nTimeIntervals; nElemICLSAC = nCCDs; - nOfInstrConstrSS = lsInstrFlag*ssInstrFlag*2*nCCDs*3; // the factor 2 comes from the AL and AC constraint equations + nOfInstrConstrSS = lsInstrFlag * ssInstrFlag * 2 * nCCDs * 3; // the factor 2 comes from the AL and AC constraint equations nElemICSS = nPixelColumns; nOfInstrConstr = nOfInstrConstrLSAL + nOfInstrConstrLSAC + nOfInstrConstrSS; - nElemIC = nOfInstrConstrLSAL*nElemICLSAL + nOfInstrConstrLSAC*nElemICLSAC + nOfInstrConstrSS*nElemICSS; - + nElemIC = nOfInstrConstrLSAL * nElemICLSAL + nOfInstrConstrLSAC * nElemICLSAC + nOfInstrConstrSS * nElemICSS; } - // Map the solved astrometric parameters, putting their indeces into an array of size nAstroPSolved // astroCoeff[0] -> parallax // astroCoeff[1] -> alpha @@ -1192,47 +1400,52 @@ int main(int argc, char **argv) { // nAstroPSolved=3 => parallax, alpha, delta, and the array is mapAstroP=[0,1,2] // nAstroPSolved=4 => alpha, delta, mu alpha, mu delta and the array is mapAstroP=[1,2,3,4] // nAstroPSolved=5 => parallax, alpha, delta, mu alpha, mu delta and the array is mapAstroP=[0,1,2,3,4] - - int LatPos=-1, LongPos=-1,MuLatPos=-1, MuLongPos=-1; - - if(nAstroPSolved) + + int LatPos = -1, LongPos = -1, MuLatPos = -1, MuLongPos = -1; + + if (nAstroPSolved) { - mapAstroP=(int *) calloc(nAstroPSolved, sizeof(int)); - switch(nAstroPSolved) { + mapAstroP = (int *)calloc(nAstroPSolved, sizeof(int)); + switch (nAstroPSolved) + { case 2: mapAstroP[0] = 1; mapAstroP[1] = 2; - LongPos=0; - LatPos=1; - MuLongPos=-1; - MuLatPos=-1; - nConstrMuLat=0; - nConstrMuLong=0; - if(extConstraint) nEqExtConstr=3; - if(barConstraint) nEqBarConstr=3; + LongPos = 0; + LatPos = 1; + MuLongPos = -1; + MuLatPos = -1; + nConstrMuLat = 0; + nConstrMuLong = 0; + if (extConstraint) + nEqExtConstr = 3; + if (barConstraint) + nEqBarConstr = 3; break; case 3: mapAstroP[0] = 0; mapAstroP[1] = 1; mapAstroP[2] = 2; - LongPos=1; - LatPos=2; - MuLongPos=-1; - MuLatPos=-1; - nConstrMuLat=0; - nConstrMuLong=0; - if(extConstraint) nEqExtConstr=3; - if(barConstraint) nEqBarConstr=3; - break; + LongPos = 1; + LatPos = 2; + MuLongPos = -1; + MuLatPos = -1; + nConstrMuLat = 0; + nConstrMuLong = 0; + if (extConstraint) + nEqExtConstr = 3; + if (barConstraint) + nEqBarConstr = 3; + break; case 4: mapAstroP[0] = 1; mapAstroP[1] = 2; mapAstroP[2] = 3; mapAstroP[3] = 4; - LongPos=0; - LatPos=1; - MuLongPos=2; - MuLatPos=3; + LongPos = 0; + LatPos = 1; + MuLongPos = 2; + MuLatPos = 3; break; case 5: mapAstroP[0] = 0; @@ -1240,861 +1453,959 @@ int main(int argc, char **argv) { mapAstroP[2] = 2; mapAstroP[3] = 3; mapAstroP[4] = 4; - LongPos=1; - LatPos=2; - MuLongPos=3; - MuLatPos=4; + LongPos = 1; + LatPos = 2; + MuLongPos = 3; + MuLatPos = 4; break; default: - if(myid==0) { + if (myid == 0) + { printf("nAstroPSolved=%d, invalid value. Aborting.\n", nAstroPSolved); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } } } // End map - /////////////////////// end broadcast section - + // Initialize the values needed to dimension the vectors nConstr = nConstrLong + nConstrLat + nConstrMuLong + nConstrMuLat; // number of constraints to be generated - - nObsxStar=nobs/nStar; - nobsOri=nobs; - if(noConstr) - { - nConstr=0; - nConstrLong=0; - nConstrLat=0; - nConstrMuLong=0; - nConstrMuLat=0; - } else - { - for(int q=0;q=nStar) - { - printf("PE=%d Error invalit Lat Constraint %ld\n",myid,constrLatId[q]); - MPI_Abort(MPI_COMM_WORLD,1); + + nObsxStar = nobs / nStar; + nobsOri = nobs; + if (noConstr) + { + nConstr = 0; + nConstrLong = 0; + nConstrLat = 0; + nConstrMuLong = 0; + nConstrMuLat = 0; + } + else + { + for (int q = 0; q < nConstrLat; q++) + if (constrLatId[q] >= nStar) + { + printf("PE=%d Error invalit Lat Constraint %ld\n", myid, constrLatId[q]); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - - for(int q=0;q=nStar) + + for (int q = 0; q < nConstrLong; q++) + if (constrLongId[q] >= nStar) { - printf("PE=%d Error invalit Long Constraint %ld\n",myid,constrLongId[q]); - MPI_Abort(MPI_COMM_WORLD,1); + printf("PE=%d Error invalit Long Constraint %ld\n", myid, constrLongId[q]); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - for(int q=0;q=nStar) + for (int q = 0; q < nConstrMuLat; q++) + if (constrMuLatId[q] >= nStar) { - printf("PE=%d Error invalit MuLat Constraint %ld\n",myid,constrMuLatId[q]); - MPI_Abort(MPI_COMM_WORLD,1); + printf("PE=%d Error invalit MuLat Constraint %ld\n", myid, constrMuLatId[q]); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - - for(int q=0;q=nStar) + + for (int q = 0; q < nConstrMuLong; q++) + if (constrMuLongId[q] >= nStar) { - printf("PE=%d Error invalit MuLat Constraint %ld\n",myid,constrMuLongId[q]); - MPI_Abort(MPI_COMM_WORLD,1); + printf("PE=%d Error invalit MuLat Constraint %ld\n", myid, constrMuLongId[q]); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - + nConstr = nConstrLong + nConstrLat + nConstrMuLong + nConstrMuLat; // number of constraints to be generated } - - - nobs+=nConstr; - nunk = ((long) nAstroParam) + nAttParam + nInstrParam + nGlobalParam; // number of unknowns (i.e. columns of the system matrix) - nparam = nAstroPSolved + nAttP + nInstrPSolved + nGlobP; // number of non-zero coefficients for each observation (i.e. for each system row) - if(nparam==0){ - printf("Abort. Empty system nparam=0 . nAstroPSolved=%d nAttP=%d nInstrPSolved=%d nGlobP=%d\n",nAstroPSolved,nAttP,nInstrPSolved,nGlobP); - MPI_Abort(MPI_COMM_WORLD,1); + + nobs += nConstr; + nunk = ((long)nAstroParam) + nAttParam + nInstrParam + nGlobalParam; // number of unknowns (i.e. columns of the system matrix) + nparam = nAstroPSolved + nAttP + nInstrPSolved + nGlobP; // number of non-zero coefficients for each observation (i.e. for each system row) + if (nparam == 0) + { + printf("Abort. Empty system nparam=0 . nAstroPSolved=%d nAttP=%d nInstrPSolved=%d nGlobP=%d\n", nAstroPSolved, nAttP, nInstrPSolved, nGlobP); + MPI_Abort(MPI_COMM_WORLD, 1); } ncoeff = nparam * nobs; // total number of non-zero coefficients of the system - - if(nobs <=nunk){ - printf("SEVERE ERROR: number of equations=%ld and number of unknown=%ld make solution unsafe\n",nobs,nunk); - MPI_Abort(MPI_COMM_WORLD,1); - exit(EXIT_FAILURE); + if (nobs <= nunk) + { + printf("SEVERE ERROR: number of equations=%ld and number of unknown=%ld make solution unsafe\n", nobs, nunk); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(EXIT_FAILURE); } - + // Start map distributed array - mapNoss=(long int *) calloc(nproc,sizeof(long int)); - mapNcoeff=(long int *) calloc(nproc,sizeof(long int)); - mapNossAfter=0; - mapNcoeffAfter=0; - mapNossBefore=0; - mapNcoeffBefore=0; - for(i=0;i=i+1) mapNoss[i]++; - mapNcoeff[i]=mapNoss[i]*nparam; - if(i= i + 1) + mapNoss[i]++; + mapNcoeff[i] = mapNoss[i] * nparam; + if (i < myid) { - mapNossBefore+=mapNoss[i]; - mapNcoeffBefore+=mapNcoeff[i]; + mapNossBefore += mapNoss[i]; + mapNcoeffBefore += mapNcoeff[i]; } - if(i>myid) + if (i > myid) { - mapNossAfter+=mapNoss[i]; - mapNcoeffAfter+=mapNcoeff[i]; + mapNossAfter += mapNoss[i]; + mapNcoeffAfter += mapNcoeff[i]; } } ////////////////// Simulating the ... of NObsxStar file - if(extConstraint){ - int * sumNObsxStar; - sumNObsxStar=(int *) calloc(nStar, sizeof(int)); - int irest=nobs % nStar; - for(int i=0;imapNossBefore && firstStarConstr==-1) firstStarConstr=i; //first star assigned in extConstr - if(counterObsxStar>=mapNossBefore+mapNoss[myid] && lastStarConstr==-1){ - lastStarConstr=i; //last star assigned in extConstr (it will be eqaul to lastrStar-1 in case of overlap) - if(counterObsxStar>(mapNossBefore+mapNoss[myid]) && myid!=(nproc-1)){ - starOverlap=1; + } + long counterObsxStar = 0; + for (int i = 0; i < nStar; i++) + { + counterObsxStar += sumNObsxStar[i]; + if (counterObsxStar > mapNossBefore && firstStarConstr == -1) + firstStarConstr = i; //first star assigned in extConstr + if (counterObsxStar >= mapNossBefore + mapNoss[myid] && lastStarConstr == -1) + { + lastStarConstr = i; //last star assigned in extConstr (it will be eqaul to lastrStar-1 in case of overlap) + if (counterObsxStar > (mapNossBefore + mapNoss[myid]) && myid != (nproc - 1)) + { + starOverlap = 1; lastStarConstr--; } break; } } - numOfExtStar=lastStarConstr-firstStarConstr+1; //number of stars computed in ext Constr - - int counterAttCols=0; - startingAttColExtConstr=0; // numero di colonna di assetto escluso l'offset: la prima è 0 per il PE0 x asse - endingAttColExtConstr=0; // numero di colonna di assetto finale l'offset: la prima è nDegFreedomAtt/nproc-1 (+1 in caso di modulo) per il PE0 x asse - int attRes=nDegFreedomAtt%nproc; - startingAttColExtConstr=(nDegFreedomAtt/nproc)*myid; - if(myidmapNossBefore && firstStarConstr==-1) firstStarConstr=i; //first star assigned in barConstr - if(counterObsxStar>=mapNossBefore+mapNoss[myid] && lastStarConstr==-1){ - lastStarConstr=i; //last star assigned in barConstr (it will be eqaul to lastrStar-1 in case of overlap) - if(counterObsxStar>(mapNossBefore+mapNoss[myid]) && myid!=(nproc-1)){ - starOverlap=1; + } + int counterObsxStar = 0; + for (int i = 0; i < nStar; i++) + { + counterObsxStar += sumNObsxStar[i]; + if (counterObsxStar > mapNossBefore && firstStarConstr == -1) + firstStarConstr = i; //first star assigned in barConstr + if (counterObsxStar >= mapNossBefore + mapNoss[myid] && lastStarConstr == -1) + { + lastStarConstr = i; //last star assigned in barConstr (it will be eqaul to lastrStar-1 in case of overlap) + if (counterObsxStar > (mapNossBefore + mapNoss[myid]) && myid != (nproc - 1)) + { + starOverlap = 1; lastStarConstr--; } break; } } - numOfBarStar=lastStarConstr-firstStarConstr+1; //number of stars computed in bar Constr - + numOfBarStar = lastStarConstr - firstStarConstr + 1; //number of stars computed in bar Constr } ////////////////////// comlsqr - - - - comlsqr.nStar=nStar; - comlsqr.nAstroP=nAstroP; - comlsqr.nAstroPSolved=nAstroPSolved; - comlsqr.nAttP=nAttP; - comlsqr.nInstrP=nInstrP; - comlsqr.nInstrPSolved=nInstrPSolved; - comlsqr.nGlobP=nGlobP; - comlsqr.mapNossBefore=mapNossBefore; - comlsqr.mapNossAfter=mapNossAfter; - comlsqr.myid=myid; - comlsqr.nproc=nproc; - comlsqr.mapNoss=mapNoss; - comlsqr.mapNcoeff=mapNcoeff; - comlsqr.multMI=multMI; - comlsqr.debugMode=debugMode; - comlsqr.noCPR=noCPR; - comlsqr.nAttParam=nAttParam; - comlsqr.extConstraint=extConstraint; - comlsqr.nEqExtConstr=nEqExtConstr; - comlsqr.numOfExtStar=numOfExtStar; - comlsqr.barConstraint=barConstraint; - comlsqr.nEqBarConstr=nEqBarConstr; - comlsqr.numOfBarStar=numOfBarStar; - comlsqr.firstStarConstr=firstStarConstr; - comlsqr.lastStarConstr=lastStarConstr; - comlsqr.numOfExtAttCol=numOfExtAttCol; - comlsqr.startingAttColExtConstr=startingAttColExtConstr; - comlsqr.setBound[0]=0; - comlsqr.setBound[1]=nAstroPSolved; - comlsqr.setBound[2]=nAstroPSolved+nAttP; - comlsqr.setBound[3]=nAstroPSolved+nAttP+nInstrPSolved; - comlsqr.nDegFreedomAtt=nDegFreedomAtt; - comlsqr.nAttParAxis=nAttParAxis; - comlsqr.nAttAxes=nAttAxes; - comlsqr.nobs=nobs; - comlsqr.lsInstrFlag=lsInstrFlag; - comlsqr.ssInstrFlag=ssInstrFlag; - comlsqr.nuInstrFlag=nuInstrFlag; - comlsqr.maInstrFlag=maInstrFlag; - comlsqr.myid=myid; - comlsqr.cCDLSAACZP=cCDLSAACZP; - comlsqr.nOfInstrConstr=nOfInstrConstr; - comlsqr.nElemIC=nElemIC; - comlsqr.nElemICLSAL=nElemICLSAL; - comlsqr.nElemICLSAC=nElemICLSAC; - comlsqr.nElemICSS=nElemICSS; - comlsqr.instrConst[0]=instrConst[0]; - comlsqr.instrConst[1]=instrConst[1]; - comlsqr.instrConst[2]=instrConst[2]; - comlsqr.instrConst[3]=instrConst[3]; - comlsqr.nInstrParam=nInstrParam; - comlsqr.nGlobalParam=nGlobalParam; + + comlsqr.nStar = nStar; + comlsqr.nAstroP = nAstroP; + comlsqr.nAstroPSolved = nAstroPSolved; + comlsqr.nAttP = nAttP; + comlsqr.nInstrP = nInstrP; + comlsqr.nInstrPSolved = nInstrPSolved; + comlsqr.nGlobP = nGlobP; + comlsqr.mapNossBefore = mapNossBefore; + comlsqr.mapNossAfter = mapNossAfter; + comlsqr.myid = myid; + comlsqr.nproc = nproc; + comlsqr.mapNoss = mapNoss; + comlsqr.mapNcoeff = mapNcoeff; + comlsqr.multMI = multMI; + comlsqr.debugMode = debugMode; + comlsqr.noCPR = noCPR; + comlsqr.nAttParam = nAttParam; + comlsqr.extConstraint = extConstraint; + comlsqr.nEqExtConstr = nEqExtConstr; + comlsqr.numOfExtStar = numOfExtStar; + comlsqr.barConstraint = barConstraint; + comlsqr.nEqBarConstr = nEqBarConstr; + comlsqr.numOfBarStar = numOfBarStar; + comlsqr.firstStarConstr = firstStarConstr; + comlsqr.lastStarConstr = lastStarConstr; + comlsqr.numOfExtAttCol = numOfExtAttCol; + comlsqr.startingAttColExtConstr = startingAttColExtConstr; + comlsqr.setBound[0] = 0; + comlsqr.setBound[1] = nAstroPSolved; + comlsqr.setBound[2] = nAstroPSolved + nAttP; + comlsqr.setBound[3] = nAstroPSolved + nAttP + nInstrPSolved; + comlsqr.nDegFreedomAtt = nDegFreedomAtt; + comlsqr.nAttParAxis = nAttParAxis; + comlsqr.nAttAxes = nAttAxes; + comlsqr.nobs = nobs; + comlsqr.lsInstrFlag = lsInstrFlag; + comlsqr.ssInstrFlag = ssInstrFlag; + comlsqr.nuInstrFlag = nuInstrFlag; + comlsqr.maInstrFlag = maInstrFlag; + comlsqr.myid = myid; + comlsqr.cCDLSAACZP = cCDLSAACZP; + comlsqr.nOfInstrConstr = nOfInstrConstr; + comlsqr.nElemIC = nElemIC; + comlsqr.nElemICLSAL = nElemICLSAL; + comlsqr.nElemICLSAC = nElemICLSAC; + comlsqr.nElemICSS = nElemICSS; + comlsqr.instrConst[0] = instrConst[0]; + comlsqr.instrConst[1] = instrConst[1]; + comlsqr.instrConst[2] = instrConst[2]; + comlsqr.instrConst[3] = instrConst[3]; + comlsqr.nInstrParam = nInstrParam; + comlsqr.nGlobalParam = nGlobalParam; + + comlsqr.nthreads= nthreads; + comlsqr.ntasks= ntasks; + ///////////////////// end buidl map distributed arrays - + // Allocate the memory for the vectors and compute the total memory allocated totmem = 0; nElements = mapNcoeff[myid]; - if (extConstraint){ - addElementextStar= (lastStarConstr-firstStarConstr+1)*nAstroPSolved; - addElementAtt=numOfExtAttCol*nAttAxes; + if (extConstraint) + { + addElementextStar = (lastStarConstr - firstStarConstr + 1) * nAstroPSolved; + addElementAtt = numOfExtAttCol * nAttAxes; } - if (barConstraint){ - addElementbarStar= (lastStarConstr-firstStarConstr+1)*nAstroPSolved; + if (barConstraint) + { + addElementbarStar = (lastStarConstr - firstStarConstr + 1) * nAstroPSolved; } - nOfElextObs=addElementextStar+addElementAtt; - nOfElBarObs=addElementbarStar; - comlsqr.nOfElextObs=nOfElextObs; - comlsqr.nOfElBarObs=nOfElBarObs; + nOfElextObs = addElementextStar + addElementAtt; + nOfElBarObs = addElementbarStar; + comlsqr.nOfElextObs = nOfElextObs; + comlsqr.nOfElBarObs = nOfElBarObs; - nElements+=nOfElextObs*nEqExtConstr+nOfElBarObs*nEqBarConstr+nElemIC; - - systemMatrix = (double *) calloc(nElements,sizeof(double)); + nElements += nOfElextObs * nEqExtConstr + nOfElBarObs * nEqBarConstr + nElemIC; + + systemMatrix = (double *)calloc(nElements, sizeof(double)); if (!systemMatrix) - exit(err_malloc("systemMatrix",myid)); - totmem += nElements*sizeof(double); - - nElements = mapNoss[myid]*multMI; - matrixIndex = (long int *) calloc(nElements, sizeof(long int)); + exit(err_malloc("systemMatrix", myid)); + totmem += nElements * sizeof(double); + + nElements = mapNoss[myid] * multMI; + matrixIndex = (long int *)calloc(nElements, sizeof(long int)); if (!matrixIndex) - exit(err_malloc("matrixIndex",myid)); - totmem += nElements* sizeof(long int); - - nElements = mapNoss[myid]*nInstrPSolved+nElemIC; - instrCol = (int *) calloc(nElements, sizeof(int)); + exit(err_malloc("matrixIndex", myid)); + totmem += nElements * sizeof(long int); + + nElements = mapNoss[myid] * nInstrPSolved + nElemIC; + instrCol = (int *)calloc(nElements, sizeof(int)); if (!instrCol) - exit(err_malloc("instrCol",myid)); - totmem += nElements* sizeof(int); + exit(err_malloc("instrCol", myid)); + totmem += nElements * sizeof(int); nElements = nOfInstrConstr; - instrConstrIlung = (int *) calloc(nElements, sizeof(int)); // it is the vectorr that for each observation (in this PE) save the INDEX for the values of instr + instrConstrIlung = (int *)calloc(nElements, sizeof(int)); // it is the vectorr that for each observation (in this PE) save the INDEX for the values of instr if (!instrConstrIlung) - exit(err_malloc("instrConstrIlung",myid)); - totmem += nElements* sizeof(int); - - + exit(err_malloc("instrConstrIlung", myid)); + totmem += nElements * sizeof(int); + nElements = mapNoss[myid]; - if(extConstraint) nElements+=nEqExtConstr; - if(barConstraint) nElements+=nEqBarConstr; - if(nOfInstrConstr>0) nElements+=nOfInstrConstr; - knownTerms = (double *) calloc(nElements, sizeof(double)); + if (extConstraint) + nElements += nEqExtConstr; + if (barConstraint) + nElements += nEqBarConstr; + if (nOfInstrConstr > 0) + nElements += nOfInstrConstr; + knownTerms = (double *)calloc(nElements, sizeof(double)); if (!knownTerms) - exit(err_malloc("knownTerms",myid)); - totmem += nElements* sizeof(double); - + exit(err_malloc("knownTerms", myid)); + totmem += nElements * sizeof(double); + ielem = 0; offsetAttParam = nAstroParam; offsetInstrParam = offsetAttParam + nAttParam; offsetGlobParam = offsetInstrParam + nInstrParam; - comlsqr.offsetAttParam=offsetAttParam; - - - long rowsxFile=nobsOri/(nproc*nfileProc); - - if(withFile) - { - char filenameSim_SM[128]="SIM_PE_SM_"; - char filenameSim_MI[128]="SIM_PE_MI_"; - char filenameSim_II[128]="SIM_PE_II_"; - char filenameSim_KT[128]="SIM_PE_KT_"; - char varpe[32]=""; - char varrows[32]=""; - sprintf(varpe,"%d",myid); - sprintf(varrows,"%ld",rowsxFile); - strcat(filenameSim_SM,varrows); - strcat(filenameSim_SM,"_"); - strcat(filenameSim_SM,varpe); - strcat(filenameSim_SM,".bin"); - strcat(filenameSim_MI,varrows); - strcat(filenameSim_MI,"_"); - strcat(filenameSim_MI,varpe); - strcat(filenameSim_MI,".bin"); - strcat(filenameSim_II,varrows); - strcat(filenameSim_II,"_"); - strcat(filenameSim_II,varpe); - strcat(filenameSim_II,".bin"); - strcat(filenameSim_KT,varrows); - strcat(filenameSim_KT,"_"); - strcat(filenameSim_KT,varpe); - strcat(filenameSim_KT,".bin"); - + comlsqr.offsetAttParam = offsetAttParam; + + long rowsxFile = nobsOri / (nproc * nfileProc); + + if (withFile) + { + char filenameSim_SM[128] = "SIM_PE_SM_"; + char filenameSim_MI[128] = "SIM_PE_MI_"; + char filenameSim_II[128] = "SIM_PE_II_"; + char filenameSim_KT[128] = "SIM_PE_KT_"; + char varpe[32] = ""; + char varrows[32] = ""; + sprintf(varpe, "%d", myid); + sprintf(varrows, "%ld", rowsxFile); + strcat(filenameSim_SM, varrows); + strcat(filenameSim_SM, "_"); + strcat(filenameSim_SM, varpe); + strcat(filenameSim_SM, ".bin"); + strcat(filenameSim_MI, varrows); + strcat(filenameSim_MI, "_"); + strcat(filenameSim_MI, varpe); + strcat(filenameSim_MI, ".bin"); + strcat(filenameSim_II, varrows); + strcat(filenameSim_II, "_"); + strcat(filenameSim_II, varpe); + strcat(filenameSim_II, ".bin"); + strcat(filenameSim_KT, varrows); + strcat(filenameSim_KT, "_"); + strcat(filenameSim_KT, varpe); + strcat(filenameSim_KT, ".bin"); + double *smsim; - smsim=(double *)calloc(rowsxFile*nparam,sizeof(double)); - fpSM=fopen(filenameSim_SM,"wb"); - fwrite(smsim,sizeof(double),rowsxFile*nparam,fpSM); + smsim = (double *)calloc(rowsxFile * nparam, sizeof(double)); + fpSM = fopen(filenameSim_SM, "wb"); + fwrite(smsim, sizeof(double), rowsxFile * nparam, fpSM); fclose(fpSM); - fpKT=fopen(filenameSim_KT,"wb"); - fwrite(smsim,sizeof(double),rowsxFile,fpKT); + fpKT = fopen(filenameSim_KT, "wb"); + fwrite(smsim, sizeof(double), rowsxFile, fpKT); fclose(fpKT); free(smsim); int *intsim; - intsim=(int *)calloc(rowsxFile*nInstrPSolved,sizeof(int)); - fpII=fopen(filenameSim_II,"wb"); - fwrite(intsim,sizeof(int),rowsxFile*nInstrPSolved,fpII); + intsim = (int *)calloc(rowsxFile * nInstrPSolved, sizeof(int)); + fpII = fopen(filenameSim_II, "wb"); + fwrite(intsim, sizeof(int), rowsxFile * nInstrPSolved, fpII); fclose(fpII); free(intsim); long *misim; - misim=(long *)calloc(rowsxFile*multMI,sizeof(long)); - fpMI=fopen(filenameSim_MI,"wb"); - fwrite(misim,sizeof(long),rowsxFile*multMI,fpMI); + misim = (long *)calloc(rowsxFile * multMI, sizeof(long)); + fpMI = fopen(filenameSim_MI, "wb"); + fwrite(misim, sizeof(long), rowsxFile * multMI, fpMI); fclose(fpMI); free(misim); long nrowToRead; - if(myid==0) printf( - "Found %d SM data files in the directory. Reading the coefficients...\n", - (nfileProc)*nproc); - timeToReadFiles=MPI_Wtime(); - - for(ii=0;ii0) + obsTotal += nObsxStar; + if (startingStar < nobsOri % nStar) + obsTotal++; + if (nConstr > 0) { - for(int q=0;q0) - { - for(int q=0;q0) - { - for(int q=0;q 0) + { + for (int q = 0; q < nConstrLat; q++) + if (constrLatId[q] == currentStar) + obsStar++; + for (int q = 0; q < nConstrLong; q++) + if (constrLongId[q] == currentStar) + obsStar++; + for (int q = 0; q < nConstrMuLat; q++) + if (constrMuLatId[q] == currentStar) + obsStar++; + for (int q = 0; q < nConstrLong; q++) + if (constrMuLongId[q] == currentStar) + obsStar++; + } + } + obsStarnow = obsStar; + if (nConstr > 0) + { + for (int q = 0; q < nConstrLat; q++) + if (constrLatId[q] == currentStar) + isConstraint++; + for (int q = 0; q < nConstrLong; q++) + if (constrLongId[q] == currentStar) + isConstraint++; + for (int q = 0; q < nConstrMuLat; q++) + if (constrMuLatId[q] == currentStar) + isConstraint++; + for (int q = 0; q < nConstrMuLong; q++) + if (constrMuLongId[q] == currentStar) + isConstraint++; + } + + int offsetConstraint = isConstraint - obsStar; // number of constraint alredy computed in the previous PE + if (offsetConstraint < 0) + offsetConstraint = 0; + + int counterStarObs = 0; + rowInFile = -1; + int changedStar = 0; + int counterConstr = 0; ///////////////////////////////// /////////// RUNNING ON ALL OBSERVATIONS ///////////////////////////////// - for(ii=0;ii=nDegFreedomAtt-nAttParAxis) lastFreedom=nDegFreedomAtt-nAttParAxis; - matrixIndex[ii*multMI+(multMI-1)]=offsetAttParam+lastFreedom; - if(lastFreedom>=endFreedom || lastFreedom>=nDegFreedomAtt-nAttParAxis) freedomReached=1; - lastFreedom+=nAttParAxis; + else + { + if (lastFreedom >= nDegFreedomAtt - nAttParAxis) + lastFreedom = nDegFreedomAtt - nAttParAxis; + matrixIndex[ii * multMI + (multMI - 1)] = offsetAttParam + lastFreedom; + if (lastFreedom >= endFreedom || lastFreedom >= nDegFreedomAtt - nAttParAxis) + freedomReached = 1; + lastFreedom += nAttParAxis; } - } else { - lastFreedom=( ( (double)rand() ) / ( ((double)RAND_MAX) ) ) * (nDegFreedomAtt-nAttParAxis+1); - if(lastFreedom>nDegFreedomAtt-nAttParAxis) lastFreedom=nDegFreedomAtt-nAttParAxis; - if((obsStar-counterStarObs)<=isConstraint) //constraint + } + else + { + lastFreedom = (((double)rand()) / (((double)RAND_MAX))) * (nDegFreedomAtt - nAttParAxis + 1); + if (lastFreedom > nDegFreedomAtt - nAttParAxis) + lastFreedom = nDegFreedomAtt - nAttParAxis; + if ((obsStar - counterStarObs) <= isConstraint) //constraint { - lastFreedom=0; - constraintFound[counterConstr][0]=currentStar/1000; - constraintFound[counterConstr][1]=rowInFile; + lastFreedom = 0; + constraintFound[counterConstr][0] = currentStar / 1000; + constraintFound[counterConstr][1] = rowInFile; counterConstr++; } - matrixIndex[ii*multMI+(multMI-1)]=offsetAttParam+lastFreedom; + matrixIndex[ii * multMI + (multMI - 1)] = offsetAttParam + lastFreedom; } ///////////// generate InstrIndex - - if(!instrFreedomReached && nInstrPSolved) - { - if((obsStar-counterStarObs)<=isConstraint) - { //constraint - for(int kk=0;kkinstrEndFreedom) instrLastFreedom=instrEndFreedom; - instrCol[ii*nInstrPSolved]=instrLastFreedom; - for(int kk=1;kk instrEndFreedom) + instrLastFreedom = instrEndFreedom; + instrCol[ii * nInstrPSolved] = instrLastFreedom; + for (int kk = 1; kk < nInstrPSolved; kk++) + instrCol[ii * nInstrPSolved + kk] = (((double)rand()) / (((double)RAND_MAX))) * (nInstrParam - 1); + if (instrLastFreedom == instrEndFreedom) + instrFreedomReached = 1; instrLastFreedom++; } - } else { - if((obsStar-counterStarObs)<=isConstraint) + } + else + { + if ((obsStar - counterStarObs) <= isConstraint) { //constraint - for(int kk=0;kkisConstraint ) - { - for(int q=0;q0) - { - if(ii!=0) offsetConstraint=0; - int foundedConstraint=(obsStar-counterStarObs)+offsetConstraint; - int itis=0; - for(int q=0;q isConstraint) + { + for (int q = 0; q < nAstroPSolved; q++) + systemMatrix[ii * nparam + q] = (((double)rand()) / RAND_MAX) * 2 - 1.0; + for (int q = 0; q < nAttP + nInstrPSolved + nGlobP; q++) + systemMatrix[ii * nparam + nAstroPSolved + q] = (((double)rand()) / RAND_MAX) * 2 - 1.0; + } + else // I add a Constraint + { + for (int q = 0; q < nAstroPSolved + nAttP + nInstrPSolved + nGlobP; q++) + systemMatrix[ii * nparam + q] = 0.; + if (nAstroPSolved > 0) + { + if (ii != 0) + offsetConstraint = 0; + int foundedConstraint = (obsStar - counterStarObs) + offsetConstraint; + int itis = 0; + for (int q = 0; q < nConstrLong; q++) + if (constrLongId[q] == currentStar) + { itis++; - if(itis==foundedConstraint) - systemMatrix[ii*nparam+LongPos]=constrLongW[q]; + if (itis == foundedConstraint) + systemMatrix[ii * nparam + LongPos] = constrLongW[q]; } - for(int q=0;q0) - { - for(int q=0;q 0) + { + for (int q = 0; q < nConstrLat; q++) + if (constrLatId[q] == currentStar) { obsStar++; isConstraint++; } - for(int q=0;q=2) randVal=0.; - if(nAstroPSolved==5 && i%nAstroPSolved==0) randVal=0.; - if(nAstroPSolved==5 && i%nAstroPSolved>2 && j<3) randVal=0.; + exit(err_malloc("attNS", myid)); + if (myid == 0) + { + for (int i = 0; i < nDegFreedomAtt; i++) + attNS[i] = (((double)rand()) / RAND_MAX) * 2 - 1.0; + } + MPI_Bcast(attNS, nDegFreedomAtt, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + for (int j = 0; j < nEqExtConstr; j++) + { + for (int i = 0; i < addElementextStar + addElementAtt; i++) + { + randVal = (((double)rand()) / RAND_MAX) * 2 - 1.0; + if (i < addElementextStar) + { + if (nAstroPSolved == 3 && i % nAstroPSolved == 0) + randVal = 0.; + if (nAstroPSolved == 4 && i % nAstroPSolved >= 2) + randVal = 0.; + if (nAstroPSolved == 5 && i % nAstroPSolved == 0) + randVal = 0.; + if (nAstroPSolved == 5 && i % nAstroPSolved > 2 && j < 3) + randVal = 0.; } - if(i>=addElementextStar){ - if(j<3) randVal=1.0; - if(j==0 || j==3){ - if(i>=addElementextStar+addElementAtt/nAttAxes) randVal=0.0; + if (i >= addElementextStar) + { + if (j < 3) + randVal = 1.0; + if (j == 0 || j == 3) + { + if (i >= addElementextStar + addElementAtt / nAttAxes) + randVal = 0.0; } - if(j==1 || j==4){ - if(i=addElementextStar+2*addElementAtt/nAttAxes) randVal=0.0; + if (j == 1 || j == 4) + { + if (i < addElementextStar + addElementAtt / nAttAxes) + randVal = 0.0; + if (i >= addElementextStar + 2 * addElementAtt / nAttAxes) + randVal = 0.0; } - if(j==2 || j==5){ - if(i=2) randVal=0.; - if(nAstroPSolved==5 && i%nAstroPSolved==0) randVal=0.; - if(nAstroPSolved==5 && i%nAstroPSolved>2 && j<3) randVal=0.; - systemMatrix[mapNcoeff[myid]+nEqExtConstr*nOfElextObs+i+j*nOfElBarObs]=randVal*barConstrW; - accumulator[j]+=randVal*barConstrW; - - } - if(!idtest) - knownTerms[mapNoss[myid]+nEqExtConstr+j]=0.; - }// j=0 - if(idtest) - MPI_Allreduce(accumulator, &knownTerms[mapNoss[myid]+nEqExtConstr], nEqBarConstr, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + accumulator = (double *)calloc(nEqBarConstr, sizeof(double)); + + for (int j = 0; j < nEqBarConstr; j++) + { + for (int i = 0; i < addElementbarStar; i++) + { + randVal = (((double)rand()) / RAND_MAX) * 2 - 1.0; + if (nAstroPSolved == 3 && i % nAstroPSolved == 0) + randVal = 0.; + if (nAstroPSolved == 4 && i % nAstroPSolved >= 2) + randVal = 0.; + if (nAstroPSolved == 5 && i % nAstroPSolved == 0) + randVal = 0.; + if (nAstroPSolved == 5 && i % nAstroPSolved > 2 && j < 3) + randVal = 0.; + systemMatrix[mapNcoeff[myid] + nEqExtConstr * nOfElextObs + i + j * nOfElBarObs] = randVal * barConstrW; + accumulator[j] += randVal * barConstrW; + } + if (!idtest) + knownTerms[mapNoss[myid] + nEqExtConstr + j] = 0.; + } // j=0 + if (idtest) + MPI_Allreduce(accumulator, &knownTerms[mapNoss[myid] + nEqExtConstr], nEqBarConstr, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); free(accumulator); } //if(barConstr - + ///////////////////////// generate instrConstr on systemMatrix // There are 1 AL + 2 AC constraint equations for each time interval for the large scale parameters (total=3*nTimeIntervals) // There are 1 AL + 1 AC constraint equations for each CCD and Legendre polinomial degree (total=2*nCCDs*3) // The equations are described by three arrays: instrCoeffConstr, instrColsConstr, instrConstrIlung, which contain the // coefficients, the column indexes and the length of the non-zero coefficients of each equation respectively. - comlsqr.offsetCMag=maInstrFlag*nCCDs; // offest=0 if maInstrFlag=0 - comlsqr.offsetCnu=comlsqr.offsetCMag+nuInstrFlag*nFoVs*nCCDs; // offest=offsetCMag if nuInstrFlag=0 - comlsqr.offsetCdelta_eta=comlsqr.offsetCnu+ssInstrFlag*nCCDs*nPixelColumns; // offest=offsetCnu if ssInstrFlag=0 - comlsqr.offsetCDelta_eta_1=comlsqr.offsetCdelta_eta+lsInstrFlag*nFoVs*nCCDs*nTimeIntervals; - comlsqr.offsetCDelta_eta_2=comlsqr.offsetCdelta_eta+lsInstrFlag*2*nFoVs*nCCDs*nTimeIntervals; - comlsqr.offsetCDelta_eta_3=comlsqr.offsetCdelta_eta+lsInstrFlag*3*nFoVs*nCCDs*nTimeIntervals; - comlsqr.offsetCdelta_zeta=comlsqr.offsetCDelta_eta_3+ssInstrFlag*nCCDs*nPixelColumns; - comlsqr.offsetCDelta_zeta_1=comlsqr.offsetCdelta_zeta+lsInstrFlag*nFoVs*nCCDs*nTimeIntervals; - comlsqr.offsetCDelta_zeta_2=comlsqr.offsetCdelta_zeta+lsInstrFlag*2*nFoVs*nCCDs*nTimeIntervals; - comlsqr.nInstrPSolved=nInstrPSolved; - - - if(myid==0 && lsInstrFlag && nElemIC!=0 ){ - double * instrCoeffConstr; - instrCoeffConstr=(double *) calloc(nElemIC, sizeof(double)); - int * instrColsConstr; - instrColsConstr=(int *) calloc(nElemIC, sizeof(int)); - - if(!computeInstrConstr(comlsqr, instrCoeffConstr, instrColsConstr, instrConstrIlung)) + comlsqr.offsetCMag = maInstrFlag * nCCDs; // offest=0 if maInstrFlag=0 + comlsqr.offsetCnu = comlsqr.offsetCMag + nuInstrFlag * nFoVs * nCCDs; // offest=offsetCMag if nuInstrFlag=0 + comlsqr.offsetCdelta_eta = comlsqr.offsetCnu + ssInstrFlag * nCCDs * nPixelColumns; // offest=offsetCnu if ssInstrFlag=0 + comlsqr.offsetCDelta_eta_1 = comlsqr.offsetCdelta_eta + lsInstrFlag * nFoVs * nCCDs * nTimeIntervals; + comlsqr.offsetCDelta_eta_2 = comlsqr.offsetCdelta_eta + lsInstrFlag * 2 * nFoVs * nCCDs * nTimeIntervals; + comlsqr.offsetCDelta_eta_3 = comlsqr.offsetCdelta_eta + lsInstrFlag * 3 * nFoVs * nCCDs * nTimeIntervals; + comlsqr.offsetCdelta_zeta = comlsqr.offsetCDelta_eta_3 + ssInstrFlag * nCCDs * nPixelColumns; + comlsqr.offsetCDelta_zeta_1 = comlsqr.offsetCdelta_zeta + lsInstrFlag * nFoVs * nCCDs * nTimeIntervals; + comlsqr.offsetCDelta_zeta_2 = comlsqr.offsetCdelta_zeta + lsInstrFlag * 2 * nFoVs * nCCDs * nTimeIntervals; + comlsqr.nInstrPSolved = nInstrPSolved; + + if (myid == 0 && lsInstrFlag && nElemIC != 0) + { + double *instrCoeffConstr; + instrCoeffConstr = (double *)calloc(nElemIC, sizeof(double)); + int *instrColsConstr; + instrColsConstr = (int *)calloc(nElemIC, sizeof(int)); + + if (!computeInstrConstr(comlsqr, instrCoeffConstr, instrColsConstr, instrConstrIlung)) { printf("SEVERE ERROR PE=0 computeInstrConstr failed\n"); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } ////////////////////////// - for(int k=0;k0) + if (seqStar <= 1 && nAstroPSolved > 0) { - printf("ERROR PE=%d Only %d star Run not allowed with this PE numbers .n",myid,seqStar); + printf("ERROR PE=%d Only %d star Run not allowed with this PE numbers .n", myid, seqStar); exit(EXIT_FAILURE); } - comlsqr.VrIdAstroPDim=seqStar; - long tempDimMax=comlsqr.VrIdAstroPDim; + comlsqr.VrIdAstroPDim = seqStar; + long tempDimMax = comlsqr.VrIdAstroPDim; long tempVrIdAstroPDimMax; - MPI_Allreduce(&tempDimMax,&tempVrIdAstroPDimMax,1,MPI_LONG,MPI_MAX,MPI_COMM_WORLD); - - comlsqr.VrIdAstroPDimMax=tempVrIdAstroPDimMax; - + MPI_Allreduce(&tempDimMax, &tempVrIdAstroPDimMax, 1, MPI_LONG, MPI_MAX, MPI_COMM_WORLD); + + comlsqr.VrIdAstroPDimMax = tempVrIdAstroPDimMax; + int **tempStarSend, **tempStarRecv; - tempStarSend=(int **) calloc(nproc,sizeof(int *)); - for(int i=0;i0) + + if (comlsqr.mapStar[myid][0] == comlsqr.mapStar[myid][1] && nAstroPSolved > 0) { - printf("PE=%d ERROR. Only one star in this PE: starting star=%d ending start=%d\n",myid,comlsqr.mapStar[myid][0],comlsqr.mapStar[myid][1]); - MPI_Abort(MPI_COMM_WORLD,1); + printf("PE=%d ERROR. Only one star in this PE: starting star=%d ending start=%d\n", myid, comlsqr.mapStar[myid][0], comlsqr.mapStar[myid][1]); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); - } - if(myid==0) - for(int i=0;id_name); + if (debugMode) + for (ii = 0; ii < nFiles; ii++) + { + printf("%s\n", namelistMI[ii]->d_name); } - if(nFiles <=0) + if (nFiles <= 0) { - printf("error on scandir n=%d\n",nFiles); - MPI_Abort(MPI_COMM_WORLD,1); + printf("error on scandir n=%d\n", nFiles); + MPI_Abort(MPI_COMM_WORLD, 1); } - + ////////////////// Writing FileConstr_GsrSolProps.dat - char fileConstr[512]=""; + char fileConstr[512] = ""; strcat(fileConstr, "FileConstr_"); strcat(fileConstr, filenameSolProps); - fpFilePosConstr=fopen(fileConstr,"w"); - for(int k=0;kd_name,constraintFound[k][1]); - } - if(debugMode) - for(int k=0;kd_name,constraintFound[k][1]); + fpFilePosConstr = fopen(fileConstr, "w"); + for (int k = 0; k < counterConstr; k++) + { + fileNum = constraintFound[k][0]; + fprintf(fpFilePosConstr, "%d %s %d\n", fileNum, namelistMI[fileNum]->d_name, constraintFound[k][1]); + } + if (debugMode) + for (int k = 0; k < counterConstr; k++) + { + fileNum = constraintFound[k][0]; + printf("PE=%d %d %s %d\n", myid, fileNum, namelistMI[fileNum]->d_name, constraintFound[k][1]); } - for(int np=1;npd_name,constraintFound[k][1]); + for (int np = 1; np < nproc; np++) + { + MPI_Recv(&counterConstr, 1, MPI_INT, np, 0, MPI_COMM_WORLD, &statusMpi); + MPI_Recv(&constraintFound[0][0], counterConstr * 2, MPI_INT, np, 1, MPI_COMM_WORLD, &statusMpi); + + for (int k = 0; k < counterConstr; k++) + { + fileNum = constraintFound[k][0]; + fprintf(fpFilePosConstr, "%d %s %d\n", fileNum, namelistMI[fileNum]->d_name, constraintFound[k][1]); } - if(debugMode) - for(int k=0;kd_name,constraintFound[k][1]); + if (debugMode) + for (int k = 0; k < counterConstr; k++) + { + fileNum = constraintFound[k][0]; + printf("PE=%d %d %s %d\n", np, fileNum, namelistMI[fileNum]->d_name, constraintFound[k][1]); } } fclose(fpFilePosConstr); - } else { + } + else + { MPI_Send(&counterConstr, 1, MPI_INT, 0, 0, MPI_COMM_WORLD); - MPI_Send(&constraintFound[0][0], counterConstr*2, MPI_INT, 0, 1, MPI_COMM_WORLD); - } - } - - VrIdAstroPDimMax=comlsqr.VrIdAstroPDimMax; - VrIdAstroPDim=comlsqr.VrIdAstroPDim; - VroffsetAttParam = VrIdAstroPDimMax*nAstroPSolved; - comlsqr.VroffsetAttParam=VroffsetAttParam; - - nunkSplit = VrIdAstroPDimMax*nAstroPSolved + nAttParam + nInstrParam + nGlobalParam; - comlsqr.nunkSplit=nunkSplit; - nElements = nunkSplit ; - preCondVect = (double *) calloc(nElements,sizeof(double)); + MPI_Send(&constraintFound[0][0], counterConstr * 2, MPI_INT, 0, 1, MPI_COMM_WORLD); + } + } + + VrIdAstroPDimMax = comlsqr.VrIdAstroPDimMax; + VrIdAstroPDim = comlsqr.VrIdAstroPDim; + VroffsetAttParam = VrIdAstroPDimMax * nAstroPSolved; + comlsqr.VroffsetAttParam = VroffsetAttParam; + + nunkSplit = VrIdAstroPDimMax * nAstroPSolved + nAttParam + nInstrParam + nGlobalParam; + comlsqr.nunkSplit = nunkSplit; + nElements = nunkSplit; + preCondVect = (double *)calloc(nElements, sizeof(double)); if (!preCondVect) - exit(err_malloc("preCondVect",myid)); - totmem += nElements*sizeof(double); - + exit(err_malloc("preCondVect", myid)); + totmem += nElements * sizeof(double); + nElements = nunkSplit; // the arrays v, w, x, and se require the same amount of memory (nunk * sizeof(double)) - vVect = (double *) calloc(nElements, sizeof(double)); + vVect = (double *)calloc(nElements, sizeof(double)); if (!vVect) - exit(err_malloc("vVect",myid)); - totmem += nElements* sizeof(double); - - wVect = (double *) calloc(nElements, sizeof(double)); + exit(err_malloc("vVect", myid)); + totmem += nElements * sizeof(double); + + + //vVect_aux_AttP = (double *)calloc(nAttParam + nInstrParam , sizeof(double)); + + + + wVect = (double *)calloc(nElements, sizeof(double)); if (!wVect) - exit(err_malloc("wVect",myid)); - totmem += nElements* sizeof(double); - - xSolution = (double *) calloc(nElements, sizeof(double)); + exit(err_malloc("wVect", myid)); + totmem += nElements * sizeof(double); + + xSolution = (double *)calloc(nElements, sizeof(double)); if (!xSolution) - exit(err_malloc("xSolution",myid)); - totmem += nElements* sizeof(double); - - standardError = (double *) calloc(nElements, sizeof(double)); + exit(err_malloc("xSolution", myid)); + totmem += nElements * sizeof(double); + + standardError = (double *)calloc(nElements, sizeof(double)); if (!standardError) - exit(err_malloc("standardError",myid)); - totmem += nElements* sizeof(double); - - totmem += nElements* sizeof(double); //dcopy+vAuxVect locally allocated on lsqr.c - - totmem=totmem/(1024*1024); // mem in MB - + exit(err_malloc("standardError", myid)); + totmem += nElements * sizeof(double); + + totmem += nElements * sizeof(double); //dcopy+vAuxVect locally allocated on lsqr.c + + totmem = totmem / (1024 * 1024); // mem in MB + // Compute and write the total memory allocated - if(myid==0) + if (myid == 0) { + printf("LOCAL %ld MB of memory allocated on each task.\n", totmem); - printf("TOTAL MB memory allocated= %ld\n", nproc*totmem ); + printf("TOTAL MB memory allocated= %ld\n", nproc * totmem); } - - + /////////////////////////////// - - - - + MPI_Barrier(MPI_COMM_WORLD); // Compute the preconditioning vector for the system columns - - - long int aIndex=0; + + long int aIndex = 0; long instrLongIndex[6]; - long iIelem=0; - - - - for (ii = 0; ii < mapNoss[myid]*multMI; ii=ii+multMI) - { - long int numOfStarPos=0; - if(nAstroPSolved>0) numOfStarPos=matrixIndex[ii]/nAstroPSolved; // number of star associated to matrixIndex[ii] - long int numOfAttPos=matrixIndex[ii+(multMI-1)]; - long int VrIdAstroPValue=-1; // - - VrIdAstroPValue=numOfStarPos-comlsqr.mapStar[myid][0]; - if(VrIdAstroPValue==-1) - { - printf("PE=%d ERROR. Can't find gsrId for precondvect.\n",myid); - MPI_Abort(MPI_COMM_WORLD,1); + long iIelem = 0; + + for (ii = 0; ii < mapNoss[myid] * multMI; ii = ii + multMI) + { + long int numOfStarPos = 0; + if (nAstroPSolved > 0) + numOfStarPos = matrixIndex[ii] / nAstroPSolved; // number of star associated to matrixIndex[ii] + long int numOfAttPos = matrixIndex[ii + (multMI - 1)]; + long int VrIdAstroPValue = -1; // + + VrIdAstroPValue = numOfStarPos - comlsqr.mapStar[myid][0]; + if (VrIdAstroPValue == -1) + { + printf("PE=%d ERROR. Can't find gsrId for precondvect.\n", myid); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - for(int ns=0;ns= nunkSplit || ncolumn < 0 ) { - printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n",myid, ncolumn, ii, + if (ncolumn >= nunkSplit || ncolumn < 0) + { + printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n", myid, ncolumn, ii, matrixIndex[ii]); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; - if(preCondVect[ncolumn]==0.0) + if (preCondVect[ncolumn] == 0.0) printf("Astrometric: preCondVect[%ld]=0.0\n", ncolumn); aIndex++; } // - for(int naxis=0;naxis= nunkSplit || ncolumn < 0 ) { - printf("ERROR. PE=%d numOfAttPos=%ld nStar*nAstroPSolved=%ld ncolumn=%ld ns=%d naxis=%d matrixIndex[ii+%d]=%ld\n",myid,numOfAttPos,nStar*nAstroPSolved,ncolumn, ns, naxis,multMI-1,matrixIndex[ii+(multMI-1)]); - MPI_Abort(MPI_COMM_WORLD,1); + ncolumn = numOfAttPos + (VroffsetAttParam - offsetAttParam) + ns + naxis * nDegFreedomAtt; + if (ncolumn >= nunkSplit || ncolumn < 0) + { + printf("ERROR. PE=%d numOfAttPos=%ld nStar*nAstroPSolved=%ld ncolumn=%ld ns=%d naxis=%d matrixIndex[ii+%d]=%ld\n", myid, numOfAttPos, nStar * nAstroPSolved, ncolumn, ns, naxis, multMI - 1, matrixIndex[ii + (multMI - 1)]); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; - if(preCondVect[ncolumn]==0.0) - printf("Attitude: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n",myid, ncolumn,aIndex,systemMatrix[aIndex]); // if aggiunto + if (preCondVect[ncolumn] == 0.0) + printf("Attitude: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n", myid, ncolumn, aIndex, systemMatrix[aIndex]); // if aggiunto aIndex++; } ///// End of Attitude preCondVect - - - if(nInstrPSolved>0) + + if (nInstrPSolved > 0) { - for(int ns=0;ns= nunkSplit || ncolumn < 0 ) + ncolumn = offsetInstrParam + (VroffsetAttParam - offsetAttParam) + instrCol[(ii / multMI) * nInstrPSolved + ns]; + if (ncolumn >= nunkSplit || ncolumn < 0) { - printf("ERROR. PE=%d ii=%ld ",myid, ii); - for(int ke=0;ke= nunkSplit || ncolumn < 0 ) { - printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii+2]=%ld\n",myid, ncolumn, ii, + ncolumn = offsetGlobParam + (VroffsetAttParam - offsetAttParam) + ns; + if (ncolumn >= nunkSplit || ncolumn < 0) + { + printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii+2]=%ld\n", myid, ncolumn, ii, matrixIndex[ii]); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; - if(preCondVect[ncolumn]==0) + if (preCondVect[ncolumn] == 0) printf("Global: preCondVect[%ld]=0.0\n", ncolumn); // if aggiunto aIndex++; } } - + ///// precondvect for extConstr - if(extConstraint){ - if(aIndex!=mapNcoeff[myid]){ - printf("PE=%d. Error on aIndex=%ld different of mapNcoeff[%d]=%ld\n",myid,aIndex,myid,mapNcoeff[myid]); + if (extConstraint) + { + if (aIndex != mapNcoeff[myid]) + { + printf("PE=%d. Error on aIndex=%ld different of mapNcoeff[%d]=%ld\n", myid, aIndex, myid, mapNcoeff[myid]); MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - for(int i=0;i0) numOfStarPos=firstStarConstr; // number of star associated to matrixIndex[ii] - long int numOfAttPos=startingAttColExtConstr; - long int VrIdAstroPValue=-1; // - - VrIdAstroPValue=numOfStarPos-comlsqr.mapStar[myid][0]; - if(VrIdAstroPValue==-1) - { - printf("PE=%d ERROR. Can't find gsrId for precondvect.\n",myid); - MPI_Abort(MPI_COMM_WORLD,1); + for (int i = 0; i < nEqExtConstr; i++) + { + long int numOfStarPos = 0; + if (nAstroPSolved > 0) + numOfStarPos = firstStarConstr; // number of star associated to matrixIndex[ii] + long int numOfAttPos = startingAttColExtConstr; + long int VrIdAstroPValue = -1; // + + VrIdAstroPValue = numOfStarPos - comlsqr.mapStar[myid][0]; + if (VrIdAstroPValue == -1) + { + printf("PE=%d ERROR. Can't find gsrId for precondvect.\n", myid); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - for(int ns=0;ns= nunkSplit || ncolumn < 0 ) { - printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n",myid, ncolumn, ii, + if (ncolumn >= nunkSplit || ncolumn < 0) + { + printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n", myid, ncolumn, ii, matrixIndex[ii]); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; - if(preCondVect[ncolumn]==0.0) + if (preCondVect[ncolumn] == 0.0) printf("Astrometric: preCondVect[%ld]=0.0\n", ncolumn); aIndex++; } // - for(int naxis=0;naxis= nunkSplit || ncolumn < 0 ) { - printf("ERROR. PE=%d ncolumn=%ld naxis=%d j=%d\n",myid,ncolumn, naxis,j); - MPI_Abort(MPI_COMM_WORLD,1); + for (int naxis = 0; naxis < nAttAxes; naxis++) + { + for (int j = 0; j < numOfExtAttCol; j++) + { + ncolumn = VrIdAstroPDimMax * nAstroPSolved + startingAttColExtConstr + j + naxis * nDegFreedomAtt; + if (ncolumn >= nunkSplit || ncolumn < 0) + { + printf("ERROR. PE=%d ncolumn=%ld naxis=%d j=%d\n", myid, ncolumn, naxis, j); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; -// printf("%d %d %d %d\n",ncolumn, aIndex,j,naxis); - if(preCondVect[ncolumn]==0.0) - printf("Attitude: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n",myid, ncolumn,aIndex,systemMatrix[aIndex]); // if aggiunto - aIndex++; + preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; + // printf("%d %d %d %d\n",ncolumn, aIndex,j,naxis); + if (preCondVect[ncolumn] == 0.0) + printf("Attitude: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n", myid, ncolumn, aIndex, systemMatrix[aIndex]); // if aggiunto + aIndex++; } } } } ///// end precondvect for extConstr ///// precondvect for barConstr - if(barConstraint){ - for(int i=0;i0) numOfStarPos=firstStarConstr; // number of star associated to matrixIndex[ii] - long int VrIdAstroPValue=-1; // - - VrIdAstroPValue=numOfStarPos-comlsqr.mapStar[myid][0]; - if(VrIdAstroPValue==-1) - { - printf("PE=%d ERROR. Can't find gsrId for precondvect.\n",myid); - MPI_Abort(MPI_COMM_WORLD,1); + if (barConstraint) + { + for (int i = 0; i < nEqBarConstr; i++) + { + long int numOfStarPos = 0; + if (nAstroPSolved > 0) + numOfStarPos = firstStarConstr; // number of star associated to matrixIndex[ii] + long int VrIdAstroPValue = -1; // + + VrIdAstroPValue = numOfStarPos - comlsqr.mapStar[myid][0]; + if (VrIdAstroPValue == -1) + { + printf("PE=%d ERROR. Can't find gsrId for precondvect.\n", myid); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - for(int ns=0;ns= nunkSplit || ncolumn < 0 ) { - printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n",myid, ncolumn, ii, + if (ncolumn >= nunkSplit || ncolumn < 0) + { + printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n", myid, ncolumn, ii, matrixIndex[ii]); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; - if(preCondVect[ncolumn]==0.0) - printf("Astrometric: preCondVect[%ld]=0.0\n", ncolumn); + if (preCondVect[ncolumn] == 0.0) + printf("Astrometric: preCondVect[%ld]=0.0\n", ncolumn); aIndex++; } // @@ -2558,291 +2895,330 @@ int main(int argc, char **argv) { ///// end precondvect for barConstr ///// precondvect for instrConstr - if(nElemIC>0){ - for(int i=0;i= nunkSplit || ncolumn < 0 ) { - printf("ERROR on instrConstr. PE=%d ncolumn=%ld i=%d\n",myid,ncolumn, i); - MPI_Abort(MPI_COMM_WORLD,1); + if (nElemIC > 0) + { + for (int i = 0; i < nElemIC; i++) + { + ncolumn = offsetInstrParam + (VroffsetAttParam - offsetAttParam) + instrCol[mapNoss[myid] * nInstrPSolved + i]; + if (ncolumn >= nunkSplit || ncolumn < 0) + { + printf("ERROR on instrConstr. PE=%d ncolumn=%ld i=%d\n", myid, ncolumn, i); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; - if(preCondVect[ncolumn]==0.0 && debugMode) - printf("Instrument: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n",myid, ncolumn,aIndex,systemMatrix[aIndex]); + if (preCondVect[ncolumn] == 0.0 && debugMode) + printf("Instrument: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n", myid, ncolumn, aIndex, systemMatrix[aIndex]); aIndex++; - } } //////////////// - - + MPI_Barrier(MPI_COMM_WORLD); - + // ACCUMULATE d // double *dcopy; - dcopy=(double *)calloc(nAttParam,sizeof(double)); - if(!dcopy) exit(err_malloc("dcopy",myid)); - mpi_allreduce(&preCondVect[ VrIdAstroPDimMax*nAstroPSolved],dcopy,(long int) nAttParam,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); - for(i=0;inAttParam && ii-VrIdAstroPDimMax*nAstroPSolvednAttParam+nInstrParam) - printf("ERROR Global ZERO column: PE=%d preCondVect[%ld]=0.0 GlobalParam=%ld \n",myid,ii, ii-(VrIdAstroPDimMax*nAstroPSolved+nAttParam+nInstrParam)); - MPI_Abort(MPI_COMM_WORLD,1); + + if (precond) + { + + if (myid == 0) + printf("Inverting preCondVect\n"); + for (ii = 0; ii < VrIdAstroPDim * nAstroPSolved; ii++) + { + if (preCondVect[ii] == 0.0) + { + if (ii - VrIdAstroPDimMax * nAstroPSolved < nAttParam) + printf("ERROR Att ZERO column: PE=%d preCondVect[%ld]=0.0 AttParam=%ld \n", myid, ii, ii - VrIdAstroPDimMax * nAstroPSolved); + if (ii - VrIdAstroPDimMax * nAstroPSolved > nAttParam && ii - VrIdAstroPDimMax * nAstroPSolved < nAttParam + nInstrParam) + printf("ERROR Instr ZERO column: PE=%d preCondVect[%ld]=0.0 InstrParam=%ld \n", myid, ii, ii - (VrIdAstroPDimMax * nAstroPSolved + nAttParam)); + if (ii - VrIdAstroPDimMax * nAstroPSolved > nAttParam + nInstrParam) + printf("ERROR Global ZERO column: PE=%d preCondVect[%ld]=0.0 GlobalParam=%ld \n", myid, ii, ii - (VrIdAstroPDimMax * nAstroPSolved + nAttParam + nInstrParam)); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ii] = 1.0 / sqrt(preCondVect[ii]); } - for (ii = VrIdAstroPDimMax*nAstroPSolved; ii < VrIdAstroPDimMax*nAstroPSolved + nAttParam + nInstrParam + nGlobalParam; ii++) { - if(preCondVect[ii]==0.0) + for (ii = VrIdAstroPDimMax * nAstroPSolved; ii < VrIdAstroPDimMax * nAstroPSolved + nAttParam + nInstrParam + nGlobalParam; ii++) + { + if (preCondVect[ii] == 0.0) { - printf("ERROR non-Astrometric ZERO column: PE=%d preCondVect[%ld]=0.0\n",myid,ii); - MPI_Abort(MPI_COMM_WORLD,1); + printf("ERROR non-Astrometric ZERO column: PE=%d preCondVect[%ld]=0.0\n", myid, ii); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ii] = 1.0 / sqrt(preCondVect[ii]); } - } else{ - if(myid==0) printf("Setting preCondVect to 1.0\n"); - for (ii = 0; ii < VrIdAstroPDim*nAstroPSolved; ii++) { - if(preCondVect[ii]==0.0) + } + else + { + if (myid == 0) + printf("Setting preCondVect to 1.0\n"); + for (ii = 0; ii < VrIdAstroPDim * nAstroPSolved; ii++) + { + if (preCondVect[ii] == 0.0) { - printf("ERROR Astrometric ZERO column: PE=%d preCondVect[%ld]=0.0 Star=%ld\n",myid,ii,ii/nAstroPSolved); - MPI_Abort(MPI_COMM_WORLD,1); + printf("ERROR Astrometric ZERO column: PE=%d preCondVect[%ld]=0.0 Star=%ld\n", myid, ii, ii / nAstroPSolved); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ii] = 1.0; } - for (ii = VrIdAstroPDimMax*nAstroPSolved; ii < VrIdAstroPDimMax*nAstroPSolved + nAttParam + nInstrParam + nGlobalParam; ii++) { - if(preCondVect[ii]==0.0) + for (ii = VrIdAstroPDimMax * nAstroPSolved; ii < VrIdAstroPDimMax * nAstroPSolved + nAttParam + nInstrParam + nGlobalParam; ii++) + { + if (preCondVect[ii] == 0.0) { - if(ii-VrIdAstroPDimMax*nAstroPSolvednAttParam && ii-VrIdAstroPDimMax*nAstroPSolvednAttParam+nInstrParam) - printf("ERROR Global ZERO column: PE=%d preCondVect[%ld]=0.0 GlobalParam=%ld \n",myid,ii, ii-(VrIdAstroPDimMax*nAstroPSolved+nAttParam+nInstrParam)); - MPI_Abort(MPI_COMM_WORLD,1); + if (ii - VrIdAstroPDimMax * nAstroPSolved < nAttParam) + printf("ERROR Att ZERO column: PE=%d preCondVect[%ld]=0.0 AttParam=%ld \n", myid, ii, ii - VrIdAstroPDimMax * nAstroPSolved); + if (ii - VrIdAstroPDimMax * nAstroPSolved > nAttParam && ii - VrIdAstroPDimMax * nAstroPSolved < nAttParam + nInstrParam) + printf("ERROR Instr ZERO column: PE=%d preCondVect[%ld]=0.0 InstrParam=%ld \n", myid, ii, ii - (VrIdAstroPDimMax * nAstroPSolved + nAttParam)); + if (ii - VrIdAstroPDimMax * nAstroPSolved > nAttParam + nInstrParam) + printf("ERROR Global ZERO column: PE=%d preCondVect[%ld]=0.0 GlobalParam=%ld \n", myid, ii, ii - (VrIdAstroPDimMax * nAstroPSolved + nAttParam + nInstrParam)); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } preCondVect[ii] = 1.0; } - } - - if(myid==0) + + if (myid == 0) { printf("Computation of the preconditioning vector finished.\n"); } - comlsqr.nvinc=0; - comlsqr.parOss=(long) nparam; - comlsqr.nunk=nunk; - comlsqr.offsetAttParam=offsetAttParam; - comlsqr.offsetInstrParam=offsetInstrParam; - comlsqr.offsetGlobParam=offsetGlobParam; - comlsqr.nAttParam=nAttParam; - comlsqr.instrConst[0]=instrConst[0]; - comlsqr.instrConst[1]=instrConst[1]; - comlsqr.instrConst[2]=instrConst[2]; - comlsqr.instrConst[3]=instrConst[3]; - comlsqr.timeCPR=timeCPR; - comlsqr.timeLimit=timeLimit; - comlsqr.itnCPR=itnCPR; - comlsqr.itnCPRstop=itnCPRstop; - comlsqr.itnLimit=itnlim; - comlsqr.Test=0; //it is not a running test but a production run - - initThread(matrixIndex,&comlsqr); - - - endTime=MPI_Wtime()-endTime; + comlsqr.nvinc = 0; + comlsqr.parOss = (long)nparam; + comlsqr.nunk = nunk; + comlsqr.offsetAttParam = offsetAttParam; + comlsqr.offsetInstrParam = offsetInstrParam; + comlsqr.offsetGlobParam = offsetGlobParam; + comlsqr.nAttParam = nAttParam; + comlsqr.instrConst[0] = instrConst[0]; + comlsqr.instrConst[1] = instrConst[1]; + comlsqr.instrConst[2] = instrConst[2]; + comlsqr.instrConst[3] = instrConst[3]; + comlsqr.timeCPR = timeCPR; + comlsqr.timeLimit = timeLimit; + comlsqr.itnCPR = itnCPR; + comlsqr.itnCPRstop = itnCPRstop; + comlsqr.itnLimit = itnlim; + comlsqr.Test = 0; //it is not a running test but a production run + + initThread(matrixIndex, &comlsqr); + + endTime = MPI_Wtime() - endTime; MPI_Barrier(MPI_COMM_WORLD); - + //////// WRITE BIN FILES and FileConstr_GsrSolProps.dat - + //////////////////////////// noiter - if(noIter){ - if(myid==0) printf("\nEnd run: -noiter option givens.\n"); - + if (noIter) + { + if (myid == 0) + printf("\nEnd run: -noiter option givens.\n"); + MPI_Finalize(); exit(EXIT_SUCCESS); - + } /////////////// MAIN CALL - + ///////////// // This function computes the product of system matrix by precondVect. This avoids to compute the produsct in aprod for each iteration. - if(myid==0 && debugMode) + if (myid == 0 && debugMode) { - printf("TEST LSQR START nobs=%ld, nunk=%ld, damp=%f, knownTerms[0]=%f, knownTerms[%ld]=%f, atol=%f btol=%f conlim=%f itnlim=%ld systemMatrix[0]=%f, systemMatrix[%ld]=%f matrixIndex[0]=%ld matrixIndex[%ld]=%ld instrCol[0]=%d instrCol[%ld]=%d preCondVect[0]=%f preCondVect[%ld]=%f preCondVect[%ld]=%f nunkSplit=%ld\n",nobs, nunk, damp, knownTerms[0], mapNoss[myid]-1,knownTerms[mapNoss[myid]-1], atol, btol,conlim,itnlim, systemMatrix[0],mapNcoeff[myid]-1,systemMatrix[mapNcoeff[myid]-1], matrixIndex[0],mapNoss[myid]*2-1, matrixIndex[mapNoss[myid]*2-1], instrCol[0], mapNoss[myid]*nInstrPSolved -1,instrCol[mapNoss[myid]*nInstrPSolved -1], preCondVect[0],VrIdAstroPDim*nAstroPSolved-1, preCondVect[VrIdAstroPDim*nAstroPSolved-1], VrIdAstroPDim*nAstroPSolved, preCondVect[VrIdAstroPDim*nAstroPSolved], nunkSplit); + printf("TEST LSQR START nobs=%ld, nunk=%ld, damp=%f, knownTerms[0]=%f, knownTerms[%ld]=%f, atol=%f btol=%f conlim=%f itnlim=%ld systemMatrix[0]=%f, systemMatrix[%ld]=%f matrixIndex[0]=%ld matrixIndex[%ld]=%ld instrCol[0]=%d instrCol[%ld]=%d preCondVect[0]=%f preCondVect[%ld]=%f preCondVect[%ld]=%f nunkSplit=%ld\n", nobs, nunk, damp, knownTerms[0], mapNoss[myid] - 1, knownTerms[mapNoss[myid] - 1], atol, btol, conlim, itnlim, systemMatrix[0], mapNcoeff[myid] - 1, systemMatrix[mapNcoeff[myid] - 1], matrixIndex[0], mapNoss[myid] * 2 - 1, matrixIndex[mapNoss[myid] * 2 - 1], instrCol[0], mapNoss[myid] * nInstrPSolved - 1, instrCol[mapNoss[myid] * nInstrPSolved - 1], preCondVect[0], VrIdAstroPDim * nAstroPSolved - 1, preCondVect[VrIdAstroPDim * nAstroPSolved - 1], VrIdAstroPDim * nAstroPSolved, preCondVect[VrIdAstroPDim * nAstroPSolved], nunkSplit); } - precondSystemMatrix(systemMatrix, preCondVect, matrixIndex,instrCol,comlsqr); + precondSystemMatrix(systemMatrix, preCondVect, matrixIndex, instrCol, comlsqr); //systemMatrix is passed to lsqr already conditioned. - if(myid==0) { + if (myid == 0) + { printf("System built, ready to call LSQR.\nPlease wait. System solution is underway..."); } //////////////////// - startTime=MPI_Wtime(); + startTime = MPI_Wtime(); seconds[1] = time(NULL); tot_sec[0] = seconds[1] - seconds[0]; - comlsqr.totSec=tot_sec[0]; - if(myid==0) printf("Starting lsqr after sec: %ld\n", tot_sec[0]); + comlsqr.totSec = tot_sec[0]; + if (myid == 0) + printf("Starting lsqr after sec: %ld\n", tot_sec[0]); chdir(wpath); - if(outputDirOpt) chdir(outputDir); // go to output dir + if (outputDirOpt) + chdir(outputDir); // go to output dir + + lsqr(nobs, nunk, damp, knownTerms, vVect, wVect, xSolution, standardError, atol, btol, conlim, - (int) itnlim, &istop, &itn, &anorm, &acond, &rnorm, &arnorm, - &xnorm, systemMatrix, matrixIndex, instrCol, instrConstrIlung,preCondVect,comlsqr); - + (int)itnlim, &istop, &itn, &anorm, &acond, &rnorm, &arnorm, + &xnorm, systemMatrix, matrixIndex, instrCol, instrConstrIlung, preCondVect, comlsqr); + /////////////////////////// MPI_Barrier(MPI_COMM_WORLD); - endTime=MPI_Wtime(); - double clockTime=MPI_Wtime(); - if(myid==0) printf("Global Time for lsqr =%f sec \n",endTime-startTime); + endTime = MPI_Wtime(); + double clockTime = MPI_Wtime(); + if (myid == 0) + printf("Global Time for lsqr =%f sec \n", endTime - startTime); seconds[2] = time(NULL); tot_sec[1] = seconds[2] - seconds[1]; - - long thetaCol=0, muthetaCol=0; + + long thetaCol = 0, muthetaCol = 0; ldiv_t res_ldiv; - long flagTheta=0, flagMuTheta=0; - - if(nAstroPSolved==2) { - thetaCol=1; - muthetaCol=0; + long flagTheta = 0, flagMuTheta = 0; + + if (nAstroPSolved == 2) + { + thetaCol = 1; + muthetaCol = 0; } - else if(nAstroPSolved==3) { - thetaCol=2; - muthetaCol=0; + else if (nAstroPSolved == 3) + { + thetaCol = 2; + muthetaCol = 0; } - else if(nAstroPSolved==4) { - thetaCol=1; - muthetaCol=3; + else if (nAstroPSolved == 4) + { + thetaCol = 1; + muthetaCol = 3; } - else if(nAstroPSolved==5) { - thetaCol=2; - muthetaCol=4; + else if (nAstroPSolved == 5) + { + thetaCol = 2; + muthetaCol = 4; } - + double epsilon; - double localSumX=0, localSumXsq=0, sumX=0, sumXsq=0, average=0, rms=0; - double dev=0, localMaxDev=0, maxDev=0; - + double localSumX = 0, localSumXsq = 0, sumX = 0, sumXsq = 0, average = 0, rms = 0; + double dev = 0, localMaxDev = 0, maxDev = 0; + ///////// Each PE runs over the its Astrometric piece - if(muthetaCol==0) - for(ii=0; iilocalMaxDev) localMaxDev=dev; + epsilon = xSolution[ii] + 1.0; + localSumX -= epsilon; + dev = fabs(epsilon); + if (dev > localMaxDev) + localMaxDev = dev; } } else { - xSolution[ii]*=preCondVect[ii]; // the corrections in theta are converted for consistency with the naming conventions in the Data Model to corrections in delta by a change of sign (Mantis Issue 0013081) - if(idtest) + xSolution[ii] *= preCondVect[ii]; // the corrections in theta are converted for consistency with the naming conventions in the Data Model to corrections in delta by a change of sign (Mantis Issue 0013081) + if (idtest) { - epsilon=xSolution[ii]-1.0; - localSumX+=epsilon; - dev=fabs(epsilon); - if(dev>localMaxDev) localMaxDev=dev; + epsilon = xSolution[ii] - 1.0; + localSumX += epsilon; + dev = fabs(epsilon); + if (dev > localMaxDev) + localMaxDev = dev; } - } - if(idtest) localSumXsq+=epsilon*epsilon; - standardError[ii]*=preCondVect[ii]; + if (idtest) + localSumXsq += epsilon * epsilon; + standardError[ii] *= preCondVect[ii]; } else - for(ii=0; iilocalMaxDev) localMaxDev=dev; + for (ii = 0; ii < VrIdAstroPDim * nAstroPSolved; ii++) + { + res_ldiv = ldiv((ii - thetaCol), nAstroPSolved); + flagTheta = res_ldiv.rem; + res_ldiv = ldiv((ii - muthetaCol), nAstroPSolved); + flagMuTheta = res_ldiv.rem; + if ((flagTheta == 0) || (flagMuTheta == 0)) + { + xSolution[ii] *= (-preCondVect[ii]); + if (idtest) + { + epsilon = xSolution[ii] + 1.0; + localSumX -= epsilon; + dev = fabs(epsilon); + if (dev > localMaxDev) + localMaxDev = dev; } - } else { - xSolution[ii]*=preCondVect[ii]; // the corrections in theta are converted for consistency with the naming conventions in the Data Model to corrections in delta by a change of sign (Mantis Issue 0013081) - if(idtest) { - epsilon=xSolution[ii]-1.0; - localSumX+=epsilon; - dev=fabs(epsilon); - if(dev>localMaxDev) localMaxDev=dev; + } + else + { + xSolution[ii] *= preCondVect[ii]; // the corrections in theta are converted for consistency with the naming conventions in the Data Model to corrections in delta by a change of sign (Mantis Issue 0013081) + if (idtest) + { + epsilon = xSolution[ii] - 1.0; + localSumX += epsilon; + dev = fabs(epsilon); + if (dev > localMaxDev) + localMaxDev = dev; } } - if(idtest) localSumXsq+=epsilon*epsilon; - standardError[ii]*=preCondVect[ii]; + if (idtest) + localSumXsq += epsilon * epsilon; + standardError[ii] *= preCondVect[ii]; } //////////// End of de-preconditioning for the Astrometric unknowns - + //////////// Then only PE=0 runs over the shared unknowns (Attitude, Instrument, and Global) - if(myid==0) - for(ii=VroffsetAttParam; iilocalMaxDev) localMaxDev=dev; - localSumXsq+=(xSolution[ii]-1.0)*(xSolution[ii]-1.0); + if (myid == 0) + for (ii = VroffsetAttParam; ii < nunkSplit; ii++) + { + xSolution[ii] *= preCondVect[ii]; + if (idtest) + { + localSumX += (xSolution[ii] - 1.0); + dev = fabs(1.0 - xSolution[ii]); + if (dev > localMaxDev) + localMaxDev = dev; + localSumXsq += (xSolution[ii] - 1.0) * (xSolution[ii] - 1.0); } - standardError[ii]*=preCondVect[ii]; + standardError[ii] *= preCondVect[ii]; } //////////// End of de-preconditioning for the shared unknowns - - + //////////////// TEST per verificare le somme di idtest.... da commentare/eliminare /* if(idtest) @@ -2854,39 +3230,72 @@ int main(int argc, char **argv) { } */ //////////////// Fine TEST - + + //print matrix + if(myid==0 ) + { + FILE * fp1; + // FILE * fp2; + + + fp1 = fopen ("./vVect.txt","w"); + // fp2 = fopen ("./matrixIndex.txt","w"); + + + for (ii = 0; ii < nunkSplit; ii++) + { + fprintf(fp1, "%16.12le \n", vVect[ii]); + // fprintf(fp2, "%16.12le \n", matrixIndex[ii]); + + + //fprintf("knownTerms[%ld]=%16.12le \n", ii,knownTerms[ii]); + } + fclose (fp1); + //fclose (fp2); + + } + + + free(systemMatrix); free(matrixIndex); free(instrCol); free(knownTerms); + + + + + + free(vVect); free(wVect); free(preCondVect); - if(idtest) + if (idtest) { - MPI_Reduce(&localSumX,&sumX,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD); - MPI_Reduce(&localSumXsq,&sumXsq,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD); - MPI_Reduce(&localMaxDev,&maxDev,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD); + MPI_Reduce(&localSumX, &sumX, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + MPI_Reduce(&localSumXsq, &sumXsq, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); + MPI_Reduce(&localMaxDev, &maxDev, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); } - - if(myid==0) + + if (myid == 0) { - if(idtest) + if (idtest) { - average=sumX/nunk; - rms=pow(sumXsq/nunk-pow(average,2.0),0.5); + average = sumX / nunk; + rms = pow(sumXsq / nunk - pow(average, 2.0), 0.5); printf("Average deviation from ID solution: %le.\n", average); printf(" rms: %le.\n", rms); printf("Maximum deviation from ID solution: %le.\n", maxDev); } } - - if(istop==1000 && wrsol==0) + + if (istop == 1000 && wrsol == 0) { - if(myid==0){ - printf("Reached limit at itn=%d. Execution stopped. CPR written!\n",itn); - ffcont=fopen("CPR_CONT","w"); + if (myid == 0) + { + printf("Reached limit at itn=%d. Execution stopped. CPR written!\n", itn); + ffcont = fopen("CPR_CONT", "w"); fclose(ffcont); } MPI_Finalize(); @@ -2894,283 +3303,304 @@ int main(int argc, char **argv) { } /////////////// ///////////// Initialize the output filename - if(istop==1000){ - sprintf(filenameAstroResults, "%s%d%s", "GsrAstroParamSolution_intermediate_",itn,".bin"); // file storing the Astrometric Parameters - sprintf(filenameAttResults, "%s%d%s","GsrAttitudeParamSolution_intermediate_",itn,".bin"); // file storing the Attitude Parameters - sprintf(filenameInstrResults, "%s%d%s","GsrInstrParamSolution_intermediate_",itn,".bin"); // file storing the Instrument Parameters - sprintf(filenameGlobalResults, "%s%d%s","GsrGlobalParamSolution_intermediate_",itn,".bin"); // file storing the Global Parameters - sprintf(filenameSolPropsFinal,"GsrFinal_intermediate_%d_%s",itn, argv[argc-1]); - } else{ - if(myid==0){ + if (istop == 1000) + { + sprintf(filenameAstroResults, "%s%d%s", "GsrAstroParamSolution_intermediate_", itn, ".bin"); // file storing the Astrometric Parameters + sprintf(filenameAttResults, "%s%d%s", "GsrAttitudeParamSolution_intermediate_", itn, ".bin"); // file storing the Attitude Parameters + sprintf(filenameInstrResults, "%s%d%s", "GsrInstrParamSolution_intermediate_", itn, ".bin"); // file storing the Instrument Parameters + sprintf(filenameGlobalResults, "%s%d%s", "GsrGlobalParamSolution_intermediate_", itn, ".bin"); // file storing the Global Parameters + sprintf(filenameSolPropsFinal, "GsrFinal_intermediate_%d_%s", itn, argv[argc - 1]); + } + else + { + if (myid == 0) + { printf("Execution finished END_FILE written!\n"); chdir(wpath); - ffcont=fopen("END_FILE","w"); + ffcont = fopen("END_FILE", "w"); fclose(ffcont); system("rm CPR_CONT"); - if(outputDirOpt) chdir(outputDir); // go to output dir - } - sprintf(filenameAstroResults, "%s", "GsrAstroParamSolution.bin"); // file storing the Astrometric Parameters - sprintf(filenameAttResults, "%s","GsrAttitudeParamSolution.bin"); // file storing the Attitude Parameters - sprintf(filenameInstrResults, "%s","GsrInstrParamSolution.bin"); // file storing the Instrument Parameters - sprintf(filenameGlobalResults, "%s","GsrGlobalParamSolution.bin"); // file storing the Global Parameters - sprintf(filenameSolPropsFinal,"GsrFinal_%s", argv[argc-1]); + if (outputDirOpt) + chdir(outputDir); // go to output dir + } + sprintf(filenameAstroResults, "%s", "GsrAstroParamSolution.bin"); // file storing the Astrometric Parameters + sprintf(filenameAttResults, "%s", "GsrAttitudeParamSolution.bin"); // file storing the Attitude Parameters + sprintf(filenameInstrResults, "%s", "GsrInstrParamSolution.bin"); // file storing the Instrument Parameters + sprintf(filenameGlobalResults, "%s", "GsrGlobalParamSolution.bin"); // file storing the Global Parameters + sprintf(filenameSolPropsFinal, "GsrFinal_%s", argv[argc - 1]); } - if(debugMode && myid==0) printf("Output %s %s %s %s %s\n",filenameSolProps,filenameAstroResults, filenameAttResults, filenameInstrResults, filenameGlobalResults); + if (debugMode && myid == 0) + printf("Output %s %s %s %s %s\n", filenameSolProps, filenameAstroResults, filenameAttResults, filenameInstrResults, filenameGlobalResults); ///////////////////////////////////////////// - + MPI_Status statusMpi; if (myid == 0) { printf("Processing finished.\n"); fflush(stdout); - for (ii = 0; ii < 10; ii++) printf("xSolution[%ld]=%16.12le \t standardError[%ld]=%16.12le\n", ii, xSolution[ii],ii,standardError[ii]); + for (ii = 0; ii < 10; ii++) + printf("xSolution[%ld]=%16.12le \t standardError[%ld]=%16.12le \n", ii, xSolution[ii], ii, standardError[ii]); } - + ////////// Writing final solution - - if(myid==0) + + if (myid == 0) { - + ////////// Writing the raw results to the files GsrFinal_GsrSolProps.dat - - - fpSolPropsFinal=fopen(filenameSolPropsFinal,"w"); + + fpSolPropsFinal = fopen(filenameSolPropsFinal, "w"); if (!fpSolPropsFinal) { - printf("Error while open %s\n",filenameSolPropsFinal); - MPI_Abort(MPI_COMM_WORLD,1); + printf("Error while open %s\n", filenameSolPropsFinal); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - fprintf(fpSolPropsFinal, "sphereId= %ld\n",sphereId) ; - fprintf(fpSolPropsFinal, "nStar= %ld\n",nStar) ; - fprintf(fpSolPropsFinal, "nAttParam= %d\n",nAttParam) ; - fprintf(fpSolPropsFinal, "nInstrParam= %d\n",nInstrParam) ; - fprintf(fpSolPropsFinal, "nInstrParamTot= %d\n",nInstrParamTot) ; - fprintf(fpSolPropsFinal, "nGlobalParam= %d\n",nGlobalParam) ; - fprintf(fpSolPropsFinal, "nAstroPSolved= %d\n",nAstroPSolved) ; - fprintf(fpSolPropsFinal, "nInstrPSolved= %d\n",nInstrPSolved) ; - fprintf(fpSolPropsFinal, "nFoVs= %d\n",nFoVs); - fprintf(fpSolPropsFinal, "nCCDs= %d\n",nCCDs); - fprintf(fpSolPropsFinal, "nPixelColumns= %d\n",nPixelColumns); - fprintf(fpSolPropsFinal, "nTimeIntervals= %d\n",nTimeIntervals); - fprintf(fpSolPropsFinal, "lSQRacond= %lf\n",acond) ; - fprintf(fpSolPropsFinal, "lSQRanorm= %lf\n",anorm) ; - fprintf(fpSolPropsFinal, "lSQRarnorm= %lf\n",arnorm) ; - fprintf(fpSolPropsFinal, "lSQRitNumb= %d\n",itn) ; - fprintf(fpSolPropsFinal, "lSQRrnorm= %lf\n",rnorm) ; - fprintf(fpSolPropsFinal, "lSQRstopCond= %d\n",istop) ; - fprintf(fpSolPropsFinal, "lSQRxnorm= %lf\n",xnorm) ; + fprintf(fpSolPropsFinal, "sphereId= %ld\n", sphereId); + fprintf(fpSolPropsFinal, "nStar= %ld\n", nStar); + fprintf(fpSolPropsFinal, "nAttParam= %d\n", nAttParam); + fprintf(fpSolPropsFinal, "nInstrParam= %d\n", nInstrParam); + fprintf(fpSolPropsFinal, "nInstrParamTot= %d\n", nInstrParamTot); + fprintf(fpSolPropsFinal, "nGlobalParam= %d\n", nGlobalParam); + fprintf(fpSolPropsFinal, "nAstroPSolved= %d\n", nAstroPSolved); + fprintf(fpSolPropsFinal, "nInstrPSolved= %d\n", nInstrPSolved); + fprintf(fpSolPropsFinal, "nFoVs= %d\n", nFoVs); + fprintf(fpSolPropsFinal, "nCCDs= %d\n", nCCDs); + fprintf(fpSolPropsFinal, "nPixelColumns= %d\n", nPixelColumns); + fprintf(fpSolPropsFinal, "nTimeIntervals= %d\n", nTimeIntervals); + fprintf(fpSolPropsFinal, "lSQRacond= %lf\n", acond); + fprintf(fpSolPropsFinal, "lSQRanorm= %lf\n", anorm); + fprintf(fpSolPropsFinal, "lSQRarnorm= %lf\n", arnorm); + fprintf(fpSolPropsFinal, "lSQRitNumb= %d\n", itn); + fprintf(fpSolPropsFinal, "lSQRrnorm= %lf\n", rnorm); + fprintf(fpSolPropsFinal, "lSQRstopCond= %d\n", istop); + fprintf(fpSolPropsFinal, "lSQRxnorm= %lf\n", xnorm); fclose(fpSolPropsFinal); - + ////////////////////// Writing GsrAstroParamSolution.bin - if(nAstroPSolved) + if (nAstroPSolved) { - long testOnStars=0; - fpAstroR=fopen(filenameAstroResults,"wb"); + long testOnStars = 0; + fpAstroR = fopen(filenameAstroResults, "wb"); if (!fpAstroR) { - printf("Error while open %s\n",filenameAstroResults); - MPI_Abort(MPI_COMM_WORLD,1); + printf("Error while open %s\n", filenameAstroResults); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - fwrite(&sphereId,sizeof(long),1,fpAstroR); - xAstro = (double *) calloc(VrIdAstroPDimMax*nAstroP, sizeof(double)); + fwrite(&sphereId, sizeof(long), 1, fpAstroR); + xAstro = (double *)calloc(VrIdAstroPDimMax * nAstroP, sizeof(double)); if (!xAstro) { printf("Error allocating xAstro\n"); - MPI_Abort(MPI_COMM_WORLD,1); + MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } - for(ii=0; ii +#include +#include +#include + +using namespace std; + +int main(int argc, char *argv[]) +{ + ifstream myfile(argv[1]); + ifstream myfile1(argv[2]); + ofstream outfile("res.txt"); + + if (myfile.is_open() && myfile1.is_open()) + { + int arrSize = 0; + double arr[1001805]; // Ideally this would be a vector, but you said array + double arr1[1001805]; // Ideally this would be a vector, but you said array + + double res[1001805]; // Ideally this would be a vector, but you said array + + + while ( true) + { + double x; + myfile >> x; + if (myfile.eof()) + break; + arr[arrSize++] = x; + } + + arrSize=0; + while (true) + { + double x; + myfile1 >> x; + if (myfile1.eof()) + break; + arr1[arrSize] = x; + + res[arrSize]=arr[arrSize]-arr1[arrSize]; + outfile<myid; +/* comlsqr->nthreads=1; #ifdef OMP comlsqr->nthreads = omp_get_max_threads(); #endif +*/ int nthreads=comlsqr->nthreads; +int ntasks=comlsqr->ntasks; + /// Prepare the structure for the division of the for cycle in aprod mode=2 -comlsqr->mapForThread=(long **) calloc(nthreads,sizeof(long *)); -for(int n=0;nmapForThread=(long **) calloc(ntasks,sizeof(long *)); +for(int n=0;nmapForThread[n]=(long *) calloc(3,sizeof(long)); -int nElements=comlsqr->mapNoss[myid]/nthreads; +int nElements=comlsqr->mapNoss[myid]/ntasks; comlsqr->mapForThread[0][0]=0; comlsqr->mapForThread[0][1]=nElements/2; comlsqr->mapForThread[0][2]=nElements; -if(comlsqr->mapNoss[myid]%nthreads>0) comlsqr->mapForThread[0][2]++; +if(comlsqr->mapNoss[myid]%ntasks>0) comlsqr->mapForThread[0][2]++; -for(int n=1;nmapForThread[n][0]=comlsqr->mapForThread[n-1][2]; comlsqr->mapForThread[n][1]=comlsqr->mapForThread[n][0]+nElements/2; comlsqr->mapForThread[n][2]=comlsqr->mapForThread[n][0]+nElements; - if(comlsqr->mapNoss[myid]%nthreads>n) comlsqr->mapForThread[n][2]++; + if(comlsqr->mapNoss[myid]%ntasks>n) comlsqr->mapForThread[n][2]++; } -comlsqr->mapForThread[nthreads-1][2]=comlsqr->mapNoss[myid]; +comlsqr->mapForThread[ntasks-1][2]=comlsqr->mapNoss[myid]; ////////////////////////////////// // Check for the NOT super-imposed stars at half cycle if(comlsqr->nAstroPSolved>0){ int smpFailed=0; - for(int n=1;nmapForThread[n-1][2]-1)]==matrixIndex[2*comlsqr->mapForThread[n][0]]) @@ -690,7 +694,7 @@ if(comlsqr->nAstroPSolved>0){ printf("UTIL: SEVERE WARNING PE=%d smpFailed\n",myid); comlsqr->mapForThread[0][0]=0; comlsqr->mapForThread[0][1]=comlsqr->mapNoss[myid]; comlsqr->mapForThread[0][2]=comlsqr->mapForThread[0][1]; - for(int n=1;nmapForThread[n][0]=comlsqr->mapNoss[myid]; comlsqr->mapForThread[n][1]=comlsqr->mapNoss[myid]; @@ -700,24 +704,28 @@ if(comlsqr->nAstroPSolved>0){ } ///// -if(comlsqr->myid==0) printf("\n\nRunning with OMP: nthreads=%d\n\n",nthreads); +if(comlsqr->myid==0) printf("\n\nRunning with OmpSs: ntasks=%d\n\n",ntasks); } +//not used void blocksAttIndep(long int *matrixIndex,struct comData *comlsqr) { //ATTENZIONE NON E? CORRETTA E NON VA MAI USATA SE nAstroPSolved=0 int myid=comlsqr->myid; - +/* comlsqr->nthreads=1; #ifdef OMP comlsqr->nthreads = omp_get_num_threads(); #endif +*/ int nthreads=comlsqr->nthreads; +int ntasks=comlsqr->ntasks; + comlsqr->nSubsetAtt=1; comlsqr->NOnSubsetAtt=1; -if(nthreads==1) +if(ntasks==1) return; int dependancyFound; @@ -729,8 +737,8 @@ dependancyFound=0; for(int i=0;inSubsetAtt;i++) //ciclo su tutti i sotto intervall { long totalEleInBlock=comlsqr->mapNoss[myid]/comlsqr->nSubsetAtt; //numero di elementi totali inclusi in tutte le thread nel blocco TBV - long totalEleInBlockThread= (comlsqr->mapNoss[myid]/comlsqr->nSubsetAtt)/nthreads; //elementi nella singola thread TBV - for(int nsubBlock=0;nsubBlockmapNoss[myid]/comlsqr->nSubsetAtt)/ntasks; //elementi nella singola thread TBV + for(int nsubBlock=0;nsubBlock nparam=nAstroPSolved+comlsqr.nAttP+comlsqr.nInstrP+comlsqr.nGlobP; nparam=nAstroPSolved+comlsqr.nAttP+comlsqr.nInstrPSolved+comlsqr.nGlobP; nLocalStar=comlsqr.mapStar[myid][1]-comlsqr.mapStar[myid][0]+1; @@ -1792,6 +1806,7 @@ struct nullSpace cknullSpace(double * systemMatrix,long * matrixIndex,double *at */ }//pragma double normLoc; + //FV_ aggiungere ntasks a cblas_dnrm2 normLoc=cblas_dnrm2(mapNoss[myid],productNS,1); double normLoc2=normLoc*normLoc; double nrmGlob; @@ -2123,3 +2138,303 @@ float simfullram(long &nStar, long &nobs, float memGlobal, int nparam, int nAttP return prevmemGB; } +double aprodM1Obs(long ix,struct comData comlsqr, double *vVect, double *systemMatrix, long int *matrixIndex, int *instrCol){ + + int myid = comlsqr.myid; + short nAstroPSolved = comlsqr.nAstroPSolved; + long localAstroMax = comlsqr.VrIdAstroPDimMax * nAstroPSolved; + long nDegFreedomAtt = comlsqr.nDegFreedomAtt; + short nAttParAxis = comlsqr.nAttParAxis; + short nAttP = comlsqr.nAttP; + short nGlobP = comlsqr.nGlobP; + short nInstrPSolved = comlsqr.nInstrPSolved; + long offsetAttParam = comlsqr.offsetAttParam; + long offsetInstrParam = comlsqr.offsetInstrParam; + long offsetGlobParam = comlsqr.offsetGlobParam; + long lset=0; + long nparam = comlsqr.parOss; + short nAttAxes = comlsqr.nAttAxes; + int multMI = comlsqr.multMI; + long miValAstro = 0; + long miValAtt = 0; + long jstartAtt = 0; + long offLocalAstro = 0; + long offLocalAtt = 0; + long offLocalInstr = 0; //Offset on Instruments + long jstartAstro = 0; + long ixInstr = 0; + int nInstrVal = 0; + double sum=0.; + offLocalInstr = offsetInstrParam + (localAstroMax - offsetAttParam); //Offset on Instruments + nInstrVal = nAstroPSolved + nAttP; + offLocalAstro = comlsqr.mapStar[myid][0] * nAstroPSolved; //Offset on my mstars + offLocalAtt = localAstroMax - offsetAttParam; //Offset on attitude + long offLocalGlob = offsetGlobParam + (localAstroMax - offsetAttParam); //Offset on GlobP + int nGlobVal = nAstroPSolved + nAttP + nInstrPSolved; + jstartAstro = miValAstro - offLocalAstro; + + ///////////////////////////////////////////////////// + /// Mode 1 Astrometric Sect + if (nAstroPSolved) + { + + lset = ix * nparam; + if (matrixIndex[multMI * ix] != miValAstro) + { + miValAstro = matrixIndex[multMI * ix]; + jstartAstro = miValAstro - offLocalAstro; + } + for (long jx = jstartAstro; jx < jstartAstro + nAstroPSolved; jx++) + { + sum = sum + systemMatrix[lset] * vVect[jx]; + lset++; + } + } + ////////////////////////////////////////////////////// + /// Mode 1 Attitude Sect + if (nAttP) + { + lset = ix * nparam + nAstroPSolved; + miValAtt = matrixIndex[multMI * ix + (multMI - 1)]; + for (int nax = 0; nax < nAttAxes; nax++) + { + jstartAtt = miValAtt + offLocalAtt + nax * nDegFreedomAtt; + for (long inpax = jstartAtt; inpax < jstartAtt + nAttParAxis; inpax++) + { + sum = sum + systemMatrix[lset] * vVect[inpax]; + lset++; + } + } + } + ////////////////////////////////////////////////////// + /// Mode 1 Instrument Sect + if (nInstrPSolved) + { + + lset = ix * nparam + nInstrVal; + long iiVal = ix * nInstrPSolved; + for (int inInstr = 0; inInstr < nInstrPSolved; inInstr++) + { + ixInstr = offLocalInstr + instrCol[iiVal + inInstr]; + sum = sum + systemMatrix[lset] * vVect[ixInstr]; + lset++; + } + } + ////////////////////////////////////////////////////// + /// Mode 1 Global sect + if (nGlobP) + { + lset = ix * nparam + nGlobVal; + for (long inGlob = offLocalGlob; inGlob < offLocalGlob + nGlobP; inGlob++) + { + sum = sum + systemMatrix[lset] * vVect[inGlob]; + lset++; + } + } + ////////////////////////////////////////////////////// + return sum; +} +void aprodM2AstroP(long ix,int nt, struct comData comlsqr, double *vVect, double *systemMatrix, long int *matrixIndex, double *knownTerms){ + + int myid = comlsqr.myid; + short nAstroPSolved = comlsqr.nAstroPSolved; + long nparam = comlsqr.parOss; + int multMI = comlsqr.multMI; + long miValAstro = 0; + long offLocalAstro = 0; + double taskLocalSum=0.; + long jstartAstro = 0; + long lset; + offLocalAstro = comlsqr.mapStar[myid][0] * nAstroPSolved; //Offset on my mstars + jstartAstro = miValAstro - offLocalAstro; + // lset=comlsqr.mapForThread[nt][0]*nparam+(ix-comlsqr.mapForThread[nt][0])*nparam; + lset=ix*nparam; + + if (matrixIndex[multMI * ix] != miValAstro) + { + miValAstro = matrixIndex[multMI * ix]; + jstartAstro = miValAstro - offLocalAstro; + } + for (long jx = jstartAstro; jx < jstartAstro + nAstroPSolved; jx++) + { + + taskLocalSum = systemMatrix[lset] * knownTerms[ix]; + //vVect[jx] = 0.000001; + if(jx == 1 && false) + { + printf("nAstroPSolved nt:%d ix:%ld vVect[jx]:%f\n\n",nt,ix,vVect[jx] ); + + } + vVect[jx] += taskLocalSum; + lset++; + } //for jx + return; +} + +void aprodM2AttP(long ix,int nt, struct comData comlsqr, double *vVect, double *systemMatrix, long int *matrixIndex, double *knownTerms){ + + int myid = comlsqr.myid; + short nAstroPSolved = comlsqr.nAstroPSolved; + long localAstroMax = comlsqr.VrIdAstroPDimMax * nAstroPSolved; + long offsetAttParam = comlsqr.offsetAttParam; + short nAttParAxis = comlsqr.nAttParAxis; + long nparam = comlsqr.parOss; + int multMI = comlsqr.multMI; + long mix; + long miValAstro = 0; + long offLocalAstro = 0; + double taskLocalSum=0.; + long jstartAstro = 0; + long lset; + long offj = (localAstroMax - offsetAttParam); + long vVix; + short nAttAxes = comlsqr.nAttAxes; + double localSum; + + long nDegFreedomAtt = comlsqr.nDegFreedomAtt; + + offLocalAstro = comlsqr.mapStar[myid][0] * nAstroPSolved; //Offset on my mstars + jstartAstro = miValAstro - offLocalAstro; + + lset = nparam * ix + nAstroPSolved; + mix = matrixIndex[multMI * ix + (multMI - 1)] + offj; + for (int ly = 0; ly < nAttAxes; ly++) + { + vVix = mix + ly * nDegFreedomAtt; + if(vVix == 1) + { + printf("nAttP nt %d %ld\n\n",nt,ix); + + } + for (int lx = 0; lx < nAttParAxis; lx++) + { + localSum = systemMatrix[lset] * knownTerms[ix]; + // #pragma omp atomic + // vVect[vVix] = 0.000001; + vVect[vVix] += localSum; + lset++; + vVix++; + } //for lx + } //for ly + return; + +} + +void aprodM2InstrP(long ix,int nt, struct comData comlsqr, double *vVect, double *systemMatrix, int *instrCol, double *knownTerms){ +int myid = comlsqr.myid; +short nAstroPSolved = comlsqr.nAstroPSolved; +long localAstroMax = comlsqr.VrIdAstroPDimMax * nAstroPSolved; +long offsetAttParam = comlsqr.offsetAttParam; +short nAttP = comlsqr.nAttP; +short nInstrPSolved = comlsqr.nInstrPSolved; +long nparam = comlsqr.parOss; +long miValAstro = 0; +long offLocalAstro = 0; +long jstartAstro = 0; +long lset; +long offsetInstrParam = comlsqr.offsetInstrParam; +long offj = offsetInstrParam + (localAstroMax - offsetAttParam); +long vVix; +long iix; + double localSum; + + +offLocalAstro = comlsqr.mapStar[myid][0] * nAstroPSolved; //Offset on my mstars +jstartAstro = miValAstro - offLocalAstro; +int offLset = nAstroPSolved + nAttP; + + lset = nparam * ix + offLset; + iix = ix * nInstrPSolved; + for (int ly = 0; ly < nInstrPSolved; ly++) + { + vVix = offj + instrCol[iix + ly]; + if(vVix == 1) + { + printf("nOfInstrConstr nt %d %ld\n\n",nt,ix); + + } + localSum = systemMatrix[lset] * knownTerms[ix]; + //#pragma omp atomic + //vVect[vVix] = 0.000001; + vVect[vVix] += localSum; + lset++; + } //for ly + return; + +} +void aprodM2GlobP(long ix,int nt, struct comData comlsqr, double *vVect, double *systemMatrix, double *knownTerms){ +int myid = comlsqr.myid; +short nAstroPSolved = comlsqr.nAstroPSolved; +long localAstroMax = comlsqr.VrIdAstroPDimMax * nAstroPSolved; +long offsetAttParam = comlsqr.offsetAttParam; +short nAttP = comlsqr.nAttP; +short nInstrPSolved = comlsqr.nInstrPSolved; +long nparam = comlsqr.parOss; +long lset; +long offsetGlobParam = comlsqr.offsetGlobParam; +short nGlobP = comlsqr.nGlobP; + +long vVix; + double localSum; + + + + int offLset = nAstroPSolved + nAttP + nInstrPSolved; + long offj = offsetGlobParam + (localAstroMax - offsetAttParam); + + lset = nparam * ix + offLset; + for (long ly = 0; ly < nGlobP; ly++) + { + vVix = offj + ly; + if(vVix == 1) + { + printf("nGlobP nt %ld\n\n",ix); + + } + localSum = systemMatrix[lset] * knownTerms[ix]; + vVect[vVix] += localSum; + lset++; + } //for ly + + + + + +} +/* + int myid = comlsqr.myid; + short nAstroPSolved = comlsqr.nAstroPSolved; + long localAstroMax = comlsqr.VrIdAstroPDimMax * nAstroPSolved; + long nDegFreedomAtt = comlsqr.nDegFreedomAtt; + short nAttParAxis = comlsqr.nAttParAxis; + short nAttP = comlsqr.nAttP; + short nGlobP = comlsqr.nGlobP; + short nInstrPSolved = comlsqr.nInstrPSolved; + long offsetAttParam = comlsqr.offsetAttParam; + long offsetInstrParam = comlsqr.offsetInstrParam; + long offsetGlobParam = comlsqr.offsetGlobParam; + long nparam = comlsqr.parOss; + short nAttAxes = comlsqr.nAttAxes; + int multMI = comlsqr.multMI; + long miValAstro = 0; + long miValAtt = 0; + long jstartAtt = 0; + long offLocalAstro = 0; + long offLocalAtt = 0; + long offLocalInstr = 0; //Offset on Instruments + double taskLocalSum=0.; + long jstartAstro = 0; + long ixInstr = 0; + int nInstrVal = 0; + double sum=0.; + long lset; + offLocalInstr = offsetInstrParam + (localAstroMax - offsetAttParam); //Offset on Instruments + nInstrVal = nAstroPSolved + nAttP; + offLocalAstro = comlsqr.mapStar[myid][0] * nAstroPSolved; //Offset on my mstars + offLocalAtt = localAstroMax - offsetAttParam; //Offset on attitude + long offLocalGlob = offsetGlobParam + (localAstroMax - offsetAttParam); //Offset on GlobP + int nGlobVal = nAstroPSolved + nAttP + nInstrPSolved; + jstartAstro = miValAstro - offLocalAstro; + lset=comlsqr.mapForThread[nt][0]*nparam+(ix-comlsqr.mapForThread[nt][0])*nparam; + + */ diff --git a/util.h b/util.h old mode 100755 new mode 100644 index d5ae47fd5d2b7632d281c4de2abd6cf427e41011..cd64880d9f75d0834692ec7646888944eedad037 --- a/util.h +++ b/util.h @@ -102,7 +102,8 @@ struct comData { int timeCPR, timeLimit, itnCPR,itnCPRstop,itnCPRend, itnLimit,itn,noCPR; long offsetCMag,offsetCnu,offsetCdelta_eta,offsetCDelta_eta_1,offsetCDelta_eta_2; long offsetCDelta_eta_3,offsetCdelta_zeta,offsetCDelta_zeta_1,offsetCDelta_zeta_2; - int nthreads; + int nthreads=1; + int ntasks; long **mapForThread; int nSubsetAtt, nSubsetInstr; int NOnSubsetAtt, NOnSubsetInstr; @@ -116,6 +117,8 @@ struct comData { char *outputDir; double nullSpaceAttfact,extConstrW,barConstrW; time_t totSec; + double **vVect_aux_AttP; + double **vVect_aux_InstP; }; @@ -209,4 +212,10 @@ double legendre(int deg, double x); int computeInstrConstr (struct comData comlsqr,double * instrCoeffConstr,int * instrColsConstr,int * instrConstrIlung); void swapInstrCoeff(double * instrCoeff, long repeat, long nrows); float simfullram(long &nStar, long &nobs, float memGlobal, int nparams, int nAttParam, int nInstrParam); +double aprodM1Obs(long ix,struct comData comlsqr, double *vVect, double *systemMatrix, long int *matrixIndex, int *instrCol); +void aprodM2AstroP(long ix,int nt, struct comData comlsqr, double *vVect, double *systemMatrix, long int *matrixIndex, double *knownTerms); +void aprodM2AttP(long ix,int nt, struct comData comlsqr, double *vVect, double *systemMatrix, long int *matrixIndex, double *knownTerms); +void aprodM2InstrP(long ix,int nt, struct comData comlsqr, double *vVect, double *systemMatrix, int *instrCol, double *knownTerms); +void aprodM2GlobP(long ix,int nt, struct comData comlsqr, double *vVect, double *systemMatrix, double *knownTerms); + #endif