Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/* lsqr.c
This C version of LSQR was first created by
Michael P Friedlander <mpf@cs.ubc.ca>
as part of his BCLS package:
http://www.cs.ubc.ca/~mpf/bcls/index.html.
The present file is maintained by
Michael Saunders <saunders@stanford.edu>
31 Aug 2007: First version of this file lsqr.c obtained from
Michael Friedlander's BCLS package, svn version number
$Revision: 273 $ $Date: 2006-09-04 15:59:04 -0700 (Mon, 04 Sep 2006) $
The stopping rules were slightly altered in that version.
They have been restored to the original rules used in the f77 LSQR.
Parallelized for ESA Gaia Mission. U. becciani A. Vecchiato 2013
*/
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <stdbool.h>
#include <math.h>
#include <time.h>
#include <mpi.h>
#include "util.h"
#include <limits.h>
//#ifdef __APPLE__
// #include <vecLib/vecLib.h>
//#else
// #include "cblas.h"
//#endif
#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);
// ---------------------------------------------------------------------
// d2norm returns sqrt( a**2 + b**2 ) with precautions
// to avoid overflow.
//
// 21 Mar 1990: First version.
// ---------------------------------------------------------------------
static double
d2norm( const double a, const double b )
{
double scale;
const double zero = 0.0;
scale = fabs( a ) + fabs( b );
if (scale == zero)
return zero;
else
return scale * sqrt( (a/scale)*(a/scale) + (b/scale)*(b/scale) );
}
static void
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;
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
//-----------------------------------------------------------------------
// 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;
// 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] =
{
{"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;
//-----------------------------------------------------------------------
// 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";
///////////// 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;
////////////////////////////////
// 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<mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr;i++){
printf("PE=%d Knowterms[%d]=%e\n",myid,i,knownTerms[i]);
}
int nStar=comlsqr.nStar;
int c1=0;
int colix=0;
for(int j=0;j<mapNoss[myid];j++)
for(int i=0; i<nAstroPSolved+nAttP;i++){
if(i==0) colix=matrixIndex[j*2];
else if(i==nAstroPSolved) colix=matrixIndex[j*2+1];
else if (i==nAstroPSolved+nAttP/nAttAxes) colix=matrixIndex[j*2+1]+nDegFreedomAtt;
else if (i==nAstroPSolved+2*nAttP/nAttAxes) colix=matrixIndex[j*2+1]+2*nDegFreedomAtt;
else colix++;
printf("PE=%d systemMatrix[%ld][%d]=%f\n",myid,mapNossBefore+j,colix,systemMatrix[c1]);
c1++;
}
for(int k=0;k<nEqExtConstr;k++)
for(int i=0; i<numOfExtStar*nAstroPSolved+numOfExtAttCol*nAttAxes;i++){
printf("PE=%d systemMatrix[%d][%ld]=%f\n",myid,nobs+k,i-mapNcoeff[myid],systemMatrix[c1]);
c1++;
}
for(int k=0;k<nEqBarConstr;k++)
for(int i=0; i<numOfBarStar*nAstroPSolved;i++){
printf("PE=%d systemMatrix[%d][%ld]=%f\n",myid,nobs+k,i-mapNcoeff[myid],systemMatrix[c1]);
c1++;
}
int counter=0;
for(int k=0;k<nOfInstrConstr;k++)
for(int j=0;j<instrConstrIlung[k];j++){
printf("PE=%d systemMatrix[%d][%ld]=%f\n",myid,nobs+nEqExtConstr+k,instrCol[mapNoss[myid]*nInstrPsolved+counter],systemMatrix[c1]);
counter++;
c1++;
}
}
if(myid==0)
{
nout=fopen(lsqrOut,"a");
if (nout != NULL)
fprintf(nout, fmt_1000,
enter, m, n, damp, wantse,
atol, conlim, btol, itnlim);
else {
printf("Error while opening %s (1)\n",lsqrOut);
MPI_Abort(MPI_COMM_WORLD,1);
exit(EXIT_FAILURE);
}
fclose(nout);
}
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;
// 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 betaLoc2=betaLoc*betaLoc;
double betaGlob;
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);
//only processor 0 accumulate and distribute the v array
#pragma omp for
for(i=0;i<localAstroMax;i++)
{
vAuxVect[i]=vVect[i];
vVect[i]=0;
}
if(myid!=0)
{
#pragma omp for
for(i=localAstroMax;i<nunkSplit;i++) vVect[i]=0;
}
aprod ( 2, m, n, vVect, knownTerms, systemMatrix, matrixIndex, instrCol,instrConstrIlung,comlsqr,&ompSec );
MPI_Barrier(MPI_COMM_WORLD);
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);
#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);
#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,(int) comlsqr.nGlobalParam,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
#pragma omp for
for(i=0;i<comlsqr.nGlobalParam;i++)
{
vVect[localAstroMax+comlsqr.nAttParam+comlsqr.nInstrParam+i]=dcopy[i];
}
free(dcopy);
if(nAstroPSolved) SumCirc(vVect,comlsqr);
#pragma omp for
for(i=0;i<localAstro;i++)
{
vVect[i]+=vAuxVect[i];
}
nAstroElements=comlsqr.mapStar[myid][1]-comlsqr.mapStar[myid][0] +1;
if(myid<nproc-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;
// 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 (myid==0)
{
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);
}
// ==================================================================
// 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);
// ==================================================================
// 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)
{
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)
{
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)
{
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);
}
}
fclose(fpCPRend);
}
if(comlsqr.itnCPRend>0 && comlsqr.itnCPRend<=itn)
{
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);
}
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<nEqExtConstr+nEqBarConstr+nOfInstrConstr;kk++){
kAuxcopy[kk]=knownTerms[mapNoss[myid]+kk];
knownTerms[mapNoss[myid]+kk]=0.;
}
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<nEqExtConstr+nEqBarConstr+nOfInstrConstr;i++)
{
knownTerms[mapNoss[myid]+i]=kcopy[i]+kAuxcopy[i];
}
free(kAuxcopy);
free(kcopy);
// beta = cblas_dnrm2 ( m, u, 1 );
if(myid==0) betaLoc=cblas_dnrm2(mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr,knownTerms,1);
else betaLoc=cblas_dnrm2(mapNoss[myid],knownTerms,1);
betaLoc2=betaLoc*betaLoc;
MPI_Allreduce(&betaLoc2, &betaGlob,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
beta=sqrt(betaGlob);
// printf("LSQR TP8 PE=%d itn=%d: alpha=%15.12lf beta=%15.12lf\n",myid,itn,alpha,beta);
// Accumulate anorm = || Bk ||
// = sqrt( sum of alpha**2 + beta**2 + damp**2 ).