V3FIT
reconstruction.f
Go to the documentation of this file.
1 !*******************************************************************************
4 !
5 ! Note separating the Doxygen comment block here so detailed decription is
6 ! found in the Module not the file.
7 !
10 !*******************************************************************************
11 
12  MODULE reconstruction
13  USE stel_kinds
14  USE guassian_process
15  USE profiler
16 
17  IMPLICIT NONE
18 
19 !*******************************************************************************
20 ! reconstruction module parameters
21 !*******************************************************************************
23  INTEGER, PARAMETER :: reconstruction_no_step_type = -1
25  INTEGER, PARAMETER :: reconstruction_sl_step_type = 0
27  INTEGER, PARAMETER :: reconstruction_lm_step_type = 1
29  INTEGER, PARAMETER :: reconstruction_seg_step_type = 2
30 
32  INTEGER, PARAMETER :: reconstruction_max_step_try = 5
33 
34 !*******************************************************************************
35 ! DERIVED-TYPE DECLARATIONS
36 ! 1) reconstruction base class
37 !
38 !*******************************************************************************
39 !-------------------------------------------------------------------------------
41 !-------------------------------------------------------------------------------
42  TYPE reconstruction_class
49  INTEGER :: step_type = reconstruction_no_step_type
52  REAL (rprec) :: step_max
54  LOGICAL :: use_central = .false.
56  INTEGER :: current_step = 0
58  INTEGER :: last_para_signal = 0
59 
62  REAL (rprec), DIMENSION(:,:), POINTER :: e => null()
65  REAL (rprec), DIMENSION(:,:), POINTER :: f => null()
67  REAL (rprec), DIMENSION(:,:), POINTER :: jacobian => null()
69  REAL (rprec), DIMENSION(:,:), POINTER :: &
70  & derived_jacobian => null()
72  REAL (rprec), DIMENSION(:,:), POINTER :: hessian => null()
74  REAL (rprec), DIMENSION(:), POINTER :: gradient => null()
76  REAL (rprec), DIMENSION(:,:), POINTER :: delta_a => null()
78  REAL (rprec), DIMENSION(:), POINTER :: j_svd_w => null()
80  REAL (rprec), DIMENSION(:,:), POINTER :: j_svd_u => null()
82  REAL (rprec), DIMENSION(:,:), POINTER :: j_svd_vt => null()
84  REAL (rprec), DIMENSION(:), POINTER :: delta_a_len => null()
86  REAL (rprec), DIMENSION(:), POINTER :: svd_w => null()
87 
89  REAL (rprec), DIMENSION(:), POINTER :: exp_g2 => null()
91  REAL (rprec), DIMENSION(:), POINTER :: step_size => null()
93  INTEGER, DIMENSION(:), POINTER :: num_sv => null()
94 
95 ! Cutoff value controls.
98  REAL (rprec) :: cut_svd
101  REAL (rprec) :: cut_eff
104  REAL (rprec) :: cut_marg_eff
107  REAL (rprec) :: cut_delta_a
110  REAL (rprec) :: cut_dg2
113  REAL (rprec) :: cut_inv_svd
114 
116  REAL (rprec), DIMENSION(:,:), POINTER :: lm_ratio => null()
117 
119  REAL (rprec), DIMENSION(:,:), POINTER :: last_values => null()
120  END TYPE
121 
122 !*******************************************************************************
123 ! INTERFACE BLOCKS
124 !*******************************************************************************
125 !-------------------------------------------------------------------------------
127 !-------------------------------------------------------------------------------
128  INTERFACE reconstruction_write_step
129  MODULE PROCEDURE reconstruction_write_step1, &
130  & reconstruction_write_step2
131  END INTERFACE
132 
133  CONTAINS
134 !*******************************************************************************
135 ! CONSTRUCTION SUBROUTINES
136 !*******************************************************************************
137 !-------------------------------------------------------------------------------
166 !-------------------------------------------------------------------------------
167  FUNCTION reconstruction_construct(num_steps, num_signals, &
168  & num_derived_parameters, &
169  & num_parameters, step_type, &
170  & step_max, cut_svd, cut_eff, &
171  & cut_marg_eff, cut_delta_a, &
172  & cut_dg2, last_para_signal, &
173  & cut_inv_svd, use_central)
174 
175  IMPLICIT NONE
176 
177 ! Declare Arguments
178  TYPE (reconstruction_class), POINTER :: reconstruction_construct
179  INTEGER, INTENT(in) :: num_steps
180  INTEGER, INTENT(in) :: num_signals
181  INTEGER, INTENT(in) :: num_derived_parameters
182  INTEGER, INTENT(in) :: num_parameters
183  CHARACTER (len=*), INTENT(in) :: step_type
184  REAL (rprec), INTENT(in) :: step_max
185  REAL (rprec), INTENT(in) :: cut_svd
186  REAL (rprec), INTENT(in) :: cut_eff
187  REAL (rprec), INTENT(in) :: cut_marg_eff
188  REAL (rprec), INTENT(in) :: cut_delta_a
189  REAL (rprec), INTENT(in) :: cut_dg2
190  INTEGER, INTENT(in) :: last_para_signal
191  REAL (rprec), INTENT(in) :: cut_inv_svd
192  LOGICAL, INTENT(in) :: use_central
193 
194 ! local variables
195  INTEGER :: num_svd
196  REAL (rprec) :: start_time
197 
198 ! Start of executable code
199  start_time = profiler_get_start_time()
200 
201  ALLOCATE(reconstruction_construct)
202 
203  num_svd = min(num_signals, num_parameters)
204 
205  SELECT CASE (trim(step_type))
206 
207  CASE ('lm')
208  reconstruction_construct%step_type = &
209  & reconstruction_lm_step_type
210  ALLOCATE(reconstruction_construct%lm_ratio(num_svd,3))
211 
212  CASE ('seg')
213  reconstruction_construct%step_type = &
214  & reconstruction_seg_step_type
215 
216  CASE DEFAULT
217  reconstruction_construct%step_type = &
218  & reconstruction_sl_step_type
219 
220  END SELECT
221 
222  reconstruction_construct%step_max = step_max
223 
224  reconstruction_construct%cut_svd = cut_svd
225  reconstruction_construct%cut_eff = cut_eff
226  reconstruction_construct%cut_marg_eff = cut_marg_eff
227  reconstruction_construct%cut_delta_a = cut_delta_a
228  reconstruction_construct%cut_dg2 = cut_dg2
229 
230  reconstruction_construct%last_para_signal = last_para_signal
231 
232  reconstruction_construct%cut_inv_svd = cut_inv_svd
233 
234  reconstruction_construct%use_central = use_central
235 
236  ALLOCATE(reconstruction_construct%e(num_signals,0:num_steps))
237  ALLOCATE(reconstruction_construct%f(num_derived_parameters, &
238  & 0:num_steps))
239  ALLOCATE(reconstruction_construct%jacobian(num_signals, &
240  & num_parameters))
241  ALLOCATE(reconstruction_construct%derived_jacobian( &
242  & num_derived_parameters,num_parameters))
243  ALLOCATE(reconstruction_construct%hessian(num_parameters, &
244  & num_parameters))
245  ALLOCATE(reconstruction_construct%gradient(num_parameters))
246  ALLOCATE(reconstruction_construct%delta_a(num_parameters, &
247  & 0:num_svd))
248  ALLOCATE(reconstruction_construct%j_svd_w(num_svd))
249  ALLOCATE(reconstruction_construct%j_svd_u(num_signals, &
250  & num_signals))
251  ALLOCATE(reconstruction_construct%j_svd_vt(num_parameters, &
252  & num_parameters))
253  ALLOCATE(reconstruction_construct%delta_a_len(0:num_svd))
254  ALLOCATE(reconstruction_construct%svd_w(num_parameters))
255  reconstruction_construct%e = 0.0
256 
257  ALLOCATE(reconstruction_construct%exp_g2(num_steps))
258  ALLOCATE(reconstruction_construct%step_size(num_steps))
259  ALLOCATE(reconstruction_construct%num_sv(num_steps))
260 
261  ALLOCATE(reconstruction_construct%last_values(4,num_signals))
262  reconstruction_construct%last_values = 0.0
263 
264  CALL profiler_set_stop_time('reconstruction_construct', &
265  & start_time)
266 
267  END FUNCTION
268 
269 !*******************************************************************************
270 ! DESTRUCTION SUBROUTINES
271 !*******************************************************************************
272 !-------------------------------------------------------------------------------
278 !-------------------------------------------------------------------------------
279  SUBROUTINE reconstruction_destruct(this)
280 
281  IMPLICIT NONE
282 
283 ! Declare Arguments
284  TYPE (reconstruction_class), POINTER :: this
285 
286 ! Start of executable code
287  this%step_type = reconstruction_no_step_type
288 
289  this%cut_svd = 0.0
290  this%cut_eff = 0.0
291  this%cut_marg_eff = 0.0
292  this%cut_delta_a = 0.0
293  this%cut_dg2 = 0.0
294 
295  this%current_step = 0
296 
297  IF (ASSOCIATED(this%e)) THEN
298  DEALLOCATE(this%e)
299  this%e => null()
300  END IF
301 
302  IF (ASSOCIATED(this%f)) THEN
303  DEALLOCATE(this%f)
304  this%f => null()
305  END IF
306 
307  IF (ASSOCIATED(this%jacobian)) THEN
308  DEALLOCATE(this%jacobian)
309  this%jacobian => null()
310  END IF
311 
312  IF (ASSOCIATED(this%derived_jacobian)) THEN
313  DEALLOCATE(this%derived_jacobian)
314  this%derived_jacobian => null()
315  END IF
316 
317  IF (ASSOCIATED(this%hessian)) THEN
318  DEALLOCATE(this%hessian)
319  this%hessian => null()
320  END IF
321 
322  IF (ASSOCIATED(this%gradient)) THEN
323  DEALLOCATE(this%gradient)
324  this%gradient => null()
325  END IF
326 
327  IF (ASSOCIATED(this%delta_a)) THEN
328  DEALLOCATE(this%delta_a)
329  this%delta_a => null()
330  END IF
331 
332  IF (ASSOCIATED(this%j_svd_w)) THEN
333  DEALLOCATE(this%j_svd_w)
334  this%j_svd_w => null()
335  END IF
336 
337  IF (ASSOCIATED(this%j_svd_u)) THEN
338  DEALLOCATE(this%j_svd_u)
339  this%j_svd_u => null()
340  END IF
341 
342  IF (ASSOCIATED(this%j_svd_vt)) THEN
343  DEALLOCATE(this%j_svd_vt)
344  this%j_svd_vt => null()
345  END IF
346 
347  IF (ASSOCIATED(this%delta_a_len)) THEN
348  DEALLOCATE(this%delta_a_len)
349  this%delta_a_len => null()
350  END IF
351 
352  IF (ASSOCIATED(this%svd_w)) THEN
353  DEALLOCATE(this%svd_w)
354  this%svd_w => null()
355  END IF
356 
357  IF (ASSOCIATED(this%exp_g2)) THEN
358  DEALLOCATE(this%exp_g2)
359  this%exp_g2 => null()
360  END IF
361 
362  IF (ASSOCIATED(this%step_size)) THEN
363  DEALLOCATE(this%step_size)
364  this%step_size => null()
365  END IF
366 
367  IF (ASSOCIATED(this%num_sv)) THEN
368  DEALLOCATE(this%num_sv)
369  this%num_sv => null()
370  END IF
371 
372  IF (ASSOCIATED(this%lm_ratio)) THEN
373  DEALLOCATE(this%lm_ratio)
374  this%lm_ratio => null()
375  END IF
376 
377  IF (ASSOCIATED(this%last_values)) THEN
378  DEALLOCATE(this%last_values)
379  this%last_values => null()
380  END IF
381 
382  DEALLOCATE(this)
383 
384  END SUBROUTINE
385 
386 !*******************************************************************************
387 ! GETTER SUBROUTINES
388 !*******************************************************************************
389 !-------------------------------------------------------------------------------
397 !-------------------------------------------------------------------------------
398  FUNCTION reconstruction_get_k_use(this, num_sv)
399 
400  IMPLICIT NONE
401 
402 ! Declare Arguments
403  INTEGER :: reconstruction_get_k_use
404  TYPE (reconstruction_class), INTENT(inout) :: this
405  INTEGER, INTENT(in) :: num_sv
406 
407 ! local variables
408  REAL (rprec), DIMENSION(:), ALLOCATABLE :: dg2exp
409  REAL (rprec), DIMENSION(:), ALLOCATABLE :: exp_eff
410  REAL (rprec), DIMENSION(:), ALLOCATABLE :: marg_exp_eff
411  INTEGER :: k
412  REAL (rprec) :: denom
413  REAL (rprec) :: largest_w
414  REAL (rprec) :: start_time
415 
416 ! Start of executable code
417  start_time = profiler_get_start_time()
418 
419  ALLOCATE(dg2exp(0:num_sv))
420  ALLOCATE(exp_eff(0:num_sv))
421  ALLOCATE(marg_exp_eff(0:num_sv))
422 
423 ! Estimate the expected changes in g^2. Equation 22 in Hanson et. al.
424 ! doi: 10.1088/0029-5515/49/7/075031
425  DO k = 0, num_sv
426  dg2exp(k) = reconstruction_get_exp_dg2(this, &
427  & this%delta_a(:,k))
428 
429  this%delta_a_len(k) = sqrt(dot_product(this%delta_a(:,k), &
430  & this%delta_a(:,k)))
431 
432 ! Equation 23 in Hanson et. al. doi: 10.1088/0029-5515/49/7/075031
433  exp_eff(k) = abs(dg2exp(k))/this%delta_a_len(k)
434  END DO
435 
436 ! Although marginal efficiencies are not strictly defined for the Steepest
437 ! Descent (index 0) and just one singular value cases, for convenience define
438 ! them here.
439  marg_exp_eff(0:1) = exp_eff(0:1)
440 
441  DO k = 2, num_sv
442  denom = this%delta_a_len(k) - this%delta_a_len(k - 1)
443  IF (denom .eq. 0.0) THEN
444  IF (this%delta_a_len(k) .eq. 0.0) THEN
445  denom = 1.0e-10_rprec
446  ELSE
447  denom = 1.0e-10_rprec*this%delta_a_len(k)
448  END IF
449  END IF
450  marg_exp_eff(k) = abs(dg2exp(k)) - abs(dg2exp(k - 1))/denom
451  END DO
452 
453 ! Check for cutoffs.
454  IF (this%j_svd_w(1) .le. 0.0) THEN
455  largest_w = 1.0
456  ELSE
457  largest_w = this%j_svd_w(1)
458  END IF
459 
460 ! Find the largest number of singular values to use. All the selection criteria
461 ! could fail so initalize it to use no sigular values.
462  reconstruction_get_k_use = 0
463  DO k = num_sv, 1, -1
464  IF ((this%j_svd_w(k)/largest_w .ge. this%cut_svd) .and. &
465  & (exp_eff(k) .ge. this%cut_eff) .and. &
466  & (marg_exp_eff(k) .ge. this%cut_marg_eff) .and. &
467  & (this%delta_a_len(k) .ge. this%cut_delta_a) .and. &
468  & (abs(dg2exp(k)) .ge. this%cut_dg2)) THEN
469  reconstruction_get_k_use = k
470  EXIT
471  END IF
472  END DO
473 
474  DEALLOCATE(dg2exp)
475  DEALLOCATE(exp_eff)
476  DEALLOCATE(marg_exp_eff)
477 
478  CALL profiler_set_stop_time('reconstruction_get_k_use', &
479  & start_time)
480 
481  END FUNCTION
482 
483 !-------------------------------------------------------------------------------
494 !-------------------------------------------------------------------------------
495  FUNCTION reconstruction_get_exp_dg2(this, delta_a)
496 
497  IMPLICIT NONE
498 
499 ! Declare Arguments
500  REAL (rprec) :: reconstruction_get_exp_dg2
501  TYPE (reconstruction_class), INTENT(in) :: this
502  REAL (rprec), DIMENSION(:), INTENT(in) :: delta_a
503 
504 ! local variables
505  REAL (rprec) :: dg2exp_lin
506  REAL (rprec) :: dg2exp_quad
507  REAL (rprec) :: start_time
508 
509 ! Start of executable code
510  start_time = profiler_get_start_time()
511 
512 ! 2e * A * da
513  dg2exp_lin = &
514  & 2.0*dot_product(this%e(:,this%current_step), &
515  & matmul(this%jacobian, delta_a))
516 
517 ! da * A^T * A * da = da * alpha * da where alpha is the hessian matrix.
518 ! Equation 14 in Hanson et. al. doi: 10.1088/0029-5515/49/7/075031
519  dg2exp_quad = dot_product(delta_a, matmul(this%hessian, delta_a))
520  reconstruction_get_exp_dg2 = dg2exp_lin - dg2exp_quad
521 
522  CALL profiler_set_stop_time('reconstruction_get_exp_dg2', &
523  & start_time)
524 
525  END FUNCTION
526 
527 !-------------------------------------------------------------------------------
536 !-------------------------------------------------------------------------------
537  FUNCTION reconstruction_get_exp_g2(this, delta_a)
538 
539  IMPLICIT NONE
540 
541 ! Declare Arguments
542  REAL (rprec) :: reconstruction_get_exp_g2
543  TYPE (reconstruction_class), INTENT(in) :: this
544  REAL (rprec), DIMENSION(:), INTENT(in) :: delta_a
545 
546 ! local variables
547  REAL (rprec) :: start_time
548 
549 ! Start of executable code
550  start_time = profiler_get_start_time()
551 
552  reconstruction_get_exp_g2 = reconstruction_get_g2(this) &
553  & - reconstruction_get_exp_dg2(this, &
554  & delta_a)
555 
556  CALL profiler_set_stop_time('reconstruction_get_exp_g2', &
557  & start_time)
558 
559  END FUNCTION
560 
561 !-------------------------------------------------------------------------------
568 !-------------------------------------------------------------------------------
569  FUNCTION reconstruction_get_g2(this)
570 
571  IMPLICIT NONE
572 
573 ! Declare Arguments
574  REAL (rprec) :: reconstruction_get_g2
575  TYPE (reconstruction_class), INTENT(in) :: this
576 
577 ! local variables
578  REAL (rprec) :: start_time
579 
580 ! Start of executable code
581  start_time = profiler_get_start_time()
582 
583  reconstruction_get_g2 = &
584  & dot_product(this%e(:,this%current_step), &
585  & this%e(:,this%current_step))
586 
587  CALL profiler_set_stop_time('reconstruction_get_g2', start_time)
588 
589  END FUNCTION
590 
591 !-------------------------------------------------------------------------------
598 !-------------------------------------------------------------------------------
599  FUNCTION reconstruction_get_lastg2(this)
600 
601  IMPLICIT NONE
602 
603 ! Declare Arguments
604  REAL (rprec) :: reconstruction_get_lastg2
605  TYPE (reconstruction_class), INTENT(in) :: this
606 
607 ! local variables
608  REAL (rprec) :: start_time
609 
610 ! Start of executable code
611  start_time = profiler_get_start_time()
612 
613  reconstruction_get_lastg2 = &
614  & dot_product(this%e(:,this%current_step - 1), &
615  & this%e(:,this%current_step - 1))
616 
617  CALL profiler_set_stop_time('reconstruction_get_lastg2', &
618  & start_time)
619  END FUNCTION
620 
621 !-------------------------------------------------------------------------------
629 !-------------------------------------------------------------------------------
630  FUNCTION reconstruction_get_dg2(this)
631 
632  IMPLICIT NONE
633 
634 ! Declare Arguments
635  REAL (rprec) :: reconstruction_get_dg2
636  TYPE (reconstruction_class), INTENT(in) :: this
637 
638 ! local variables
639  REAL (rprec) :: start_time
640 
641 ! Start of executable code
642  start_time = profiler_get_start_time()
643 
644  reconstruction_get_dg2 = reconstruction_get_lastg2(this) &
645  & - reconstruction_get_g2(this)
646 
647  CALL profiler_set_stop_time('reconstruction_get_dg2', start_time)
648  END FUNCTION
649 
650 !*******************************************************************************
651 ! UTILITY SUBROUTINES
652 !*******************************************************************************
653 !-------------------------------------------------------------------------------
667 !-------------------------------------------------------------------------------
668  FUNCTION reconstruction_eval_e(this, signals, a_model, gaussp, &
669  & eq_steps, iou, eq_comm)
670 
671  IMPLICIT NONE
672 
673 ! Declare Arguments
674  LOGICAL :: reconstruction_eval_e
675  TYPE (reconstruction_class), INTENT(inout) :: this
676  TYPE (signal_pointer), DIMENSION(:), INTENT(inout) :: signals
677  TYPE (model_class), POINTER :: a_model
678  TYPE (gaussp_class_pointer), DIMENSION(:), INTENT(inout) :: &
679  & gaussp
680  INTEGER, INTENT(inout) :: eq_steps
681  INTEGER, INTENT(in) :: iou
682  INTEGER, INTENT(in) :: eq_comm
683 
684 ! local variables
685  INTEGER :: i
686  REAL (rprec) :: start_time
687 
688 ! Start of executable code
689  start_time = profiler_get_start_time()
690 
691 ! Converge the equilibrium. this%e is initialized to zero in the constructor.
692  IF (model_converge(a_model, eq_steps, iou, eq_comm, 'All')) THEN
693  DO i = 1, SIZE(gaussp)
694  CALL gaussp_set_profile(gaussp(i)%p, a_model)
695  END DO
696 
697 !$OMP PARALLEL DO
698 !$OMP& SCHEDULE(DYNAMIC)
699 !$OMP& DEFAULT(SHARED)
700 !$OMP& PRIVATE(i)
701  DO i = 1, this%last_para_signal
702  this%e(i,this%current_step) = &
703  & signals(i)%p%get_e(a_model, .false., &
704  & this%last_values(:,i))
705  END DO
706 !$OMP END PARALLEL DO
707 
708  DO i = this%last_para_signal + 1, SIZE(signals)
709  this%e(i,this%current_step) = &
710  & signals(i)%p%get_e(a_model, .false., &
711  & this%last_values(:,i))
712  END DO
713 
714 ! All signals are computed reset the model state.
715  a_model%state_flags = model_state_all_off
716 
717  ! Write out the step information.
718  reconstruction_eval_e = .true.
719  ELSE
720  reconstruction_eval_e = .false.
721  END IF
722 
723  CALL profiler_set_stop_time('reconstruction_eval_e', start_time)
724 
725  END FUNCTION
726 
727 !-------------------------------------------------------------------------------
736 !-------------------------------------------------------------------------------
737  SUBROUTINE reconstruction_eval_f(this, derived_params, a_model)
738 
739  IMPLICIT NONE
740 
741 ! Declare Arguments
742  TYPE (reconstruction_class), INTENT(inout) :: this
743  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: derived_params
744  TYPE (model_class), INTENT(inout) :: a_model
745 
746 ! local variables
747  INTEGER :: i
748  REAL (rprec) :: start_time
749 
750 ! Start of executable code
751  start_time = profiler_get_start_time()
752 
753  DO i = 1, SIZE(derived_params)
754  this%f(i,this%current_step) = &
755  & param_get_value(derived_params(i)%p, a_model)
756  END DO
757 
758  CALL profiler_set_stop_time('reconstruction_eval_f', start_time)
759 
760  END SUBROUTINE
761 
762 !-------------------------------------------------------------------------------
779 !-------------------------------------------------------------------------------
780  SUBROUTINE reconstruction_eval_jacobians(this, signals, &
781  & derived_params, locks, &
782  & a_model, gaussp, params, &
783  & eq_steps, iou, &
784  & recon_comm, eq_comm)
785  USE xstuff
786 
787  IMPLICIT NONE
788 
789 ! Declare Arguments
790  TYPE (reconstruction_class), INTENT(inout) :: this
791  TYPE (signal_pointer), DIMENSION(:), INTENT(inout) :: signals
792  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: &
793  & derived_params
794  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: locks
795  TYPE (model_class), POINTER :: a_model
796  TYPE (gaussp_class_pointer), DIMENSION(:), INTENT(inout) :: &
797  & gaussp
798  TYPE (param_pointer), DIMENSION(:), INTENT(inout) :: params
799  INTEGER, INTENT(inout) :: eq_steps
800  INTEGER, INTENT(in) :: iou
801  INTEGER, INTENT(in) :: recon_comm
802  INTEGER, INTENT(in) :: eq_comm
803 
804 ! local variables
805  INTEGER :: i, j
806  REAL (rprec) :: param_value
807  REAL (rprec) :: start_time
808  INTEGER :: mpi_size
809  INTEGER :: mpi_rank
810 #if defined(MPI_OPT)
811  INTEGER :: error
812 #endif
813  INTEGER, DIMENSION(:), ALLOCATABLE :: stride
814  INTEGER, DIMENSION(:), ALLOCATABLE :: offset
815  INTEGER :: minsize
816  INTEGER :: extra
817 
818 ! Start of
819  start_time = profiler_get_start_time()
820 
821 #if defined(MPI_OPT)
822 ! Set the params first since they can change the state of the converged
823 ! variable. Syning the model afterwards will bring the internal state up to
824 ! sync.
825  DO i = 1, SIZE(params)
826  CALL param_sync_value(params(i)%p, a_model, recon_comm, &
827  & eq_comm, this%use_central)
828  END DO
829  CALL model_sync_state(a_model, recon_comm)
830  CALL reconstruction_sync_state(this, recon_comm)
831 
832 ! Set the locked values.
833  DO i = 1, SIZE(locks)
834  CALL param_set_lock_value(locks(i)%p, a_model, eq_comm)
835  END DO
836 
837  CALL mpi_comm_rank(recon_comm, mpi_rank, error)
838  CALL mpi_comm_size(recon_comm, mpi_size, error)
839 
840 #else
841  mpi_rank = 0
842  mpi_size = 1
843 #endif
844 
845  ALLOCATE(stride(mpi_size))
846  ALLOCATE(offset(mpi_size))
847 
848 ! Determine the offsets and strides that each uses. This determined the array
849 ! indices that each process computes the jacobian for. The number of parameters
850 ! might not be evenly dividable by the number of process. As a result. need to
851 ! determine the minimum size and add extra work for some of the processes.
852  minsize = SIZE(params)/mpi_size
853  extra = mod(SIZE(params), mpi_size)
854 
855  stride = minsize
856  stride(1:extra) = stride(1:extra) + 1
857  offset = 0
858  DO i = 2, mpi_size
859  offset(i) = offset(i - 1) + stride(i - 1)
860  END DO
861 
862  CALL model_save_state(a_model)
863 
864 ! Loop over the parameters, but only do iterations for the rank as determined
865 ! by the stride and the offset. Strides and offsets are determined using C
866 ! array notation since the MPI libraries use C conventions. Make sure the
867 ! offset is incremented by one for the initial index.
868  DO j = offset(mpi_rank + 1) + 1, offset(mpi_rank + 1) + &
869  & stride(mpi_rank + 1)
870 
871 ! Save the initial parameter value and the equilirbium state.
872  param_value = param_get_value(params(j)%p, a_model)
873 
874 ! Increment the parameter. param_increment checks the bounds and adjust the
875 ! step size to keep the parameter in range.
876  CALL param_increment(params(j)%p, a_model, eq_comm, &
877  & this%use_central)
878 ! Set the locked values.
879  DO i = 1, SIZE(locks)
880  CALL param_set_lock_value(locks(i)%p, a_model, eq_comm)
881  END DO
882 
883 ! Check that model was able to converge. If the model converged proceed as
884 ! normal otherwise set the jth index of the jacobian to zero to avoid a step
885 ! in that parameter direction.
886  IF (model_converge(a_model, eq_steps, iou, eq_comm, &
887  & param_get_name(params(j)%p, a_model))) THEN
888  DO i = 1, SIZE(gaussp)
889  CALL gaussp_set_profile(gaussp(i)%p, a_model)
890  END DO
891 
892 ! Compute the normalized jacobian. d e_i/d a_j . Equation 12 in Hanson et. al.
893 ! doi:10.1088/0029-5515/49/7/075031
894 !
895 ! a_j = p_j/|pi_j| (1)
896 !
897 ! Where p is the jth parameter, and pi is the jth parameter vrnc. For the
898 ! finite difference, the partial derivative becomes
899 !
900 ! d e_i/d a_j = (e'_i - e_i)/(a'_j - a_j) (2)
901 !
902 ! Where the primed components are values for the incremented value.
903 !
904 ! a'_j = (p_j + pi_j)/|pi_j| (3)
905 !
906 ! Subsituting equation 2 into 3, reduces to
907 !
908 ! d e_i/d a_j = (e'_i - e_i)*|pi_j|/pi_j (4)
909 !
910 ! In this coding e_i = Sqrt(signal_get_g2)
911 ! and |pi_j|/pi_j = SIGN(1,delta_param)
912 !
913 ! For the derived parameter jacobian. d f_i/d a_j . This will parallel the
914 ! signal jacobian with the exception that the parameters aren't normalized.
915 ! Using equation (3) above the derived jacobian becomes
916 !
917 ! d f_i/d a_j = (f'_i - f_i)*|pi_j|/pi_j (5)
918 !
919 ! Again in this coding |pi_j|/pi_j = SIGN(1,delta_param)
920 
921  IF (this%use_central) THEN
922 ! Certian signals cannot be parallelized. Compute all the signals up to that
923 ! point in parallel then continue in serial. Store half the finite difference
924 ! in the jacobian. Central difference only takes steps in a single direction so
925 ! there is no need to normaize the direction.
926 !$OMP PARALLEL DO
927 !$OMP& SCHEDULE(DYNAMIC)
928 !$OMP& DEFAULT(SHARED)
929 !$OMP& PRIVATE(i)
930  DO i = 1, this%last_para_signal
931  this%jacobian(i,j) = &
932  & signals(i)%p%get_e(a_model, .false., &
933  & this%last_values(:,i))
934  END DO
935 !$OMP END PARALLEL DO
936  DO i = this%last_para_signal + 1, SIZE(signals)
937  this%jacobian(i,j) = &
938  & signals(i)%p%get_e(a_model, .false., &
939  & this%last_values(:,i))
940  END DO
941 
942 ! Derived Jacobian.
943  DO i = 1, SIZE(derived_params)
944  this%derived_jacobian(i,j) = &
945  & param_get_value(derived_params(i)%p, a_model)
946  END DO
947 
948 ! Reset the parameters and equilibrium to prepare for the backward step.
949  CALL param_set_value(params(j)%p, a_model, param_value, &
950  & eq_comm, this%use_central)
951  CALL model_reset_state(a_model)
952 ! Set the locked values.
953  DO i = 1, SIZE(locks)
954  CALL param_set_lock_value(locks(i)%p, a_model, &
955  & eq_comm)
956  END DO
957 
958 ! Decrement the parameter.
959  CALL param_decrement(params(j)%p, a_model, eq_comm)
960 
961  IF (model_converge(a_model, eq_steps, iou, eq_comm, &
962  & param_get_name(params(j)%p, &
963  & a_model))) THEN
964  DO i = 1, SIZE(gaussp)
965  CALL gaussp_set_profile(gaussp(i)%p, a_model)
966  END DO
967 
968 ! Find other half of the finite difference. Negate the jacobian to obtain a
969 ! direction down hill.
970 !$OMP PARALLEL DO
971 !$OMP& SCHEDULE(DYNAMIC)
972 !$OMP& DEFAULT(SHARED)
973 !$OMP& PRIVATE(i)
974  DO i = 1, this%last_para_signal
975  this%jacobian(i,j) = &
976  & -this%jacobian(i,j) + &
977  & signals(i)%p%get_e(a_model, .false., &
978  & this%last_values(:,i))
979  END DO
980 !$OMP END PARALLEL DO
981  DO i = this%last_para_signal + 1, SIZE(signals)
982  this%jacobian(i,j) = &
983  & -this%jacobian(i,j) + &
984  & signals(i)%p%get_e(a_model, .false., &
985  & this%last_values(:,i))
986  END DO
987 
988 ! Derived Jacobian.
989  DO i = 1, SIZE(derived_params)
990  this%derived_jacobian(i,j) = &
991  & -this%derived_jacobian(i,j) + &
992  & param_get_value(derived_params(i)%p, a_model)
993  END DO
994 
995  ELSE
996 ! Warn the user that the model failed to converge on the jth parameter.
997  WRITE (*,1000) j
998 
999  this%jacobian(:,j) = 0.0
1000  this%derived_jacobian(:,j) = 0.0
1001  END IF
1002  ELSE
1003 
1004 ! Certian signals cannot be parallelized. Compute all the signals up to that
1005 ! point in parallel then continue in serial. When using a single sided
1006 ! difference, the step may be taken in either direction. Normalize by the sign
1007 ! of the parameter delta.
1008 !$OMP PARALLEL DO
1009 !$OMP& SCHEDULE(DYNAMIC)
1010 !$OMP& DEFAULT(SHARED)
1011 !$OMP& PRIVATE(i)
1012  DO i = 1, this%last_para_signal
1013  this%jacobian(i,j) = &
1014  & -(signals(i)%p%get_e(a_model, .false., &
1015  & this%last_values(:,i)) - &
1016  & this%e(i,this%current_step))
1017  END DO
1018 !$OMP END PARALLEL DO
1019  DO i = this%last_para_signal + 1, SIZE(signals)
1020  this%jacobian(i,j) = &
1021  & -(signals(i)%p%get_e(a_model, .false., &
1022  & this%last_values(:,i)) - &
1023  & this%e(i,this%current_step))
1024  END DO
1025  this%jacobian(:,j) = this%jacobian(:,j) &
1026  & * sign(1.0_rprec, &
1027  & params(j)%p%recon%delta)
1028 
1029 ! Derived Jacobian.
1030  DO i = 1, SIZE(derived_params)
1031  this%derived_jacobian(i,j) = &
1032  & -(param_get_value(derived_params(i)%p, a_model) &
1033  & - this%f(i,this%current_step))
1034  END DO
1035  this%derived_jacobian(:,j) = &
1036  & this%derived_jacobian(:,j) * &
1037  & sign(1.0_rprec, params(j)%p%recon%delta)
1038  END IF
1039  ELSE
1040 ! Warn the user that the model failed to converge on the jth parameter.
1041  WRITE (*,1000) j
1042 
1043  this%jacobian(:,j) = 0.0
1044  this%derived_jacobian(:,j) = 0.0
1045  END IF
1046 
1047 ! Restore the parameter and the equilibrium to it's initial value. Need to sync
1048 ! the value of param delta back to the parent process.
1049  CALL param_set_value(params(j)%p, a_model, param_value, &
1050  & eq_comm, this%use_central)
1051  CALL model_reset_state(a_model)
1052 ! Set the locked values.
1053  DO i = 1, SIZE(locks)
1054  CALL param_set_lock_value(locks(i)%p, a_model, eq_comm)
1055  END DO
1056  END DO
1057 ! At this point, the computed jacobian rows should be returned to the parent
1058 ! process.
1059 #if defined(MPI_OPT)
1060  CALL reconstruction_sync_parent(this, params, SIZE(signals), &
1061  & SIZE(derived_params), stride, &
1062  & offset, recon_comm)
1063 #endif
1064 
1065  DEALLOCATE(stride)
1066  DEALLOCATE(offset)
1067 
1068  CALL profiler_set_stop_time('reconstruction_eval_jacobians', &
1069  & start_time)
1070 
1071 1000 FORMAT('Warning: Model failed to converge in Jacobian for ', &
1072  & 'index ',i3,'.')
1073 
1074  END SUBROUTINE
1075 
1076 !-------------------------------------------------------------------------------
1093 !-------------------------------------------------------------------------------
1094  FUNCTION reconstruction_eval_step(this, signals, derived_params, &
1095  & locks, a_model, gaussp, params, &
1096  & eq_steps, iou, recon_comm, &
1097  & eq_comm)
1098 
1099  IMPLICIT NONE
1100 
1101 ! Declare Arguments
1102  INTEGER :: reconstruction_eval_step
1103  TYPE (reconstruction_class), INTENT(inout) :: this
1104  TYPE (signal_pointer), DIMENSION(:), INTENT(inout) :: signals
1105  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: &
1106  & derived_params
1107  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: locks
1108  TYPE (model_class), POINTER :: a_model
1109  TYPE (gaussp_class_pointer), DIMENSION(:), INTENT(inout) :: &
1110  & gaussp
1111  TYPE (param_pointer), DIMENSION(:), INTENT(inout) :: params
1112  INTEGER, INTENT(inout) :: eq_steps
1113  INTEGER, INTENT(in) :: iou
1114  INTEGER, INTENT(in) :: recon_comm
1115  INTEGER, INTENT(in) :: eq_comm
1116 
1117 ! local variables
1118  INTEGER :: i
1119  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: temp_jacobian
1120  REAL (rprec), DIMENSION(:), ALLOCATABLE :: temp_work
1121  REAL (rprec), DIMENSION(:), ALLOCATABLE :: temp_work2
1122  REAL (rprec), DIMENSION(:), ALLOCATABLE :: temp_work3
1123  INTEGER :: svd_status
1124  REAL (rprec) :: start_time
1125 #if defined(MPI_OPT)
1126  INTEGER :: error
1127 #endif
1128 
1129 ! Start of executable code
1130  start_time = profiler_get_start_time()
1131 
1132 ! Tell the child processes to begin the jacobian task.
1133 #if defined(MPI_OPT)
1134  CALL mpi_bcast(mpi_jacobian_task, 1, mpi_integer, 0, recon_comm, &
1135  & error)
1136  CALL mpi_bcast(eq_steps, 1, mpi_integer, 0, recon_comm, error)
1137 #endif
1138 
1139  CALL reconstruction_eval_jacobians(this, signals, derived_params, &
1140  & locks, a_model, gaussp, params, &
1141  & eq_steps, iou, recon_comm, &
1142  & eq_comm)
1143 
1144 ! Compute the positive definite Hessian, matrix. Equation 14 in Hanson et. al.
1145 ! doi:10.1088/0029-5515/49/7/075031
1146 !
1147 ! alpha = A^T * A (1)
1148 !
1149 ! Numerical Recipes Third edition Equation 15.5.7 . In here, the second
1150 ! derivatives are assumed to be small and ignored.
1151  this%hessian = matmul(transpose(this%jacobian), this%jacobian)
1152 
1153 ! Compute the gradient vector in parameter space. Equation 13 in Hanson et. al.
1154 ! doi:10.1088/0029-5515/49/7/075031
1155 !
1156 ! beta = A^T * e (2)
1157  this%gradient = matmul(transpose(this%jacobian), &
1158  & this%e(:,this%current_step))
1159 
1160 ! Compute the steepest descent step size in normalized parameter space.
1161 !
1162 ! da = beta * beta * beta / (beta * alpha * beta) (3)
1163  this%delta_a(:,0) = this%gradient*dot_product(this%gradient, &
1164  & this%gradient) &
1165  & / dot_product(this%gradient, &
1166  & matmul(this%hessian, &
1167  & this%gradient))
1168 
1169 ! Singular value decomposition of the jacobian. Need to allocate a new jacobian
1170 ! matrix because dgesvd will destroy the one passed to it. dgesvd is from
1171 ! LAPACK
1172  ALLOCATE(temp_jacobian(SIZE(signals),SIZE(params)))
1173  temp_jacobian = this%jacobian
1174 
1175  ALLOCATE(temp_work(5*max(SIZE(signals), SIZE(params))))
1176  temp_work = 0.0
1177 
1178  CALL dgesvd('All', 'All', SIZE(signals), SIZE(params), &
1179  & temp_jacobian, SIZE(signals), this%j_svd_w, &
1180  & this%j_svd_u, SIZE(signals), this%j_svd_vt, &
1181  & SIZE(params), temp_work, SIZE(temp_work), svd_status)
1182  CALL assert_eq(0, svd_status, 'reconstruction_eval_step: ' // &
1183  & 'dgesvd problem')
1184 
1185  DEALLOCATE(temp_work)
1186  DEALLOCATE(temp_jacobian)
1187 
1188 ! Reuse the temp work array.
1189  ALLOCATE(temp_work(SIZE(this%j_svd_w)))
1190  ALLOCATE(temp_work2(SIZE(signals)))
1191  ALLOCATE(temp_work3(SIZE(params)))
1192 
1193 ! Define the inverse singular values.
1194  WHERE (this%j_svd_w .gt. 0.0)
1195  temp_work = 1.0/this%j_svd_w
1196  ELSE WHERE
1197  temp_work = 0.0
1198  END WHERE
1199 
1200 ! U^T * e
1201  temp_work2 = matmul(transpose(this%j_svd_u), &
1202  & this%e(:,this%current_step))
1203 
1204 ! Perform pseudo-inversion for successive numbers of singular values retained.
1205  temp_work3 = 0.0
1206  DO i = 1, SIZE(this%j_svd_w)
1207  temp_work3(i) = temp_work(i)*temp_work2(i)
1208  this%delta_a(:,i) = matmul(transpose(this%j_svd_vt), &
1209  & temp_work3)
1210  END DO
1211 
1212  DEALLOCATE(temp_work)
1213  DEALLOCATE(temp_work2)
1214  DEALLOCATE(temp_work3)
1215 
1216  reconstruction_eval_step = &
1217  & reconstruction_get_k_use(this, SIZE(this%j_svd_w))
1218 
1219  CALL profiler_set_stop_time('reconstruction_eval_step', &
1220  & start_time)
1221 
1222  END FUNCTION
1223 
1224 !-------------------------------------------------------------------------------
1234 !-------------------------------------------------------------------------------
1235  FUNCTION reconstruction_sl_step(this, k_use, step_size, delta_a)
1236 
1237  IMPLICIT NONE
1238 
1239 ! Declare Arguments
1240  REAL (rprec) :: reconstruction_sl_step
1241  TYPE (reconstruction_class), INTENT(inout) :: this
1242  INTEGER, INTENT(in) :: k_use
1243  REAL (rprec), INTENT(in) :: step_size
1244  REAL (rprec), DIMENSION(:), INTENT(out) :: delta_a
1245 
1246 ! local variables
1247  REAL (rprec) :: start_time
1248 
1249 ! Start of executable code
1250  start_time = profiler_get_start_time()
1251 
1252  IF ((step_size .ge. 0.0) .and. &
1253  & (this%delta_a_len(k_use) .gt. step_size)) THEN
1254  delta_a = this%step_max/this%delta_a_len(k_use) &
1255  & * this%delta_a(:,k_use)
1256  reconstruction_sl_step = sqrt(dot_product(delta_a, delta_a))
1257  ELSE
1258  delta_a = this%delta_a(:,k_use)
1259  reconstruction_sl_step = this%delta_a_len(k_use)
1260  END IF
1261 
1262  CALL profiler_set_stop_time('reconstruction_sl_step', start_time)
1263 
1264  END FUNCTION
1265 
1266 !-------------------------------------------------------------------------------
1276 !-------------------------------------------------------------------------------
1277  FUNCTION reconstruction_seg_step(this, k_use, step_size, delta_a)
1278 
1279  IMPLICIT NONE
1280 
1281 ! Declare Arguments
1282  REAL (rprec) :: reconstruction_seg_step
1283  TYPE (reconstruction_class), INTENT(inout) :: this
1284  INTEGER, INTENT(in) :: k_use
1285  REAL (rprec), INTENT(in) :: step_size
1286  REAL (rprec), DIMENSION(:), INTENT(out) :: delta_a
1287 
1288 ! local variables
1289  INTEGER :: k, k_min
1290  REAL (rprec) :: ysq, zsq, ydz
1291  REAL (rprec) :: aa, bb, cc, xx
1292  REAL (rprec) :: discriminant
1293  REAL (rprec) :: start_time
1294 
1295 ! Start of executable code
1296  start_time = profiler_get_start_time()
1297 
1298  IF ((step_size .ge. 0.0) .and. &
1299  & (this%delta_a_len(k_use) .gt. step_size)) THEN
1300 ! First, determine k_min - the largest k value where the step is smaller than
1301 ! step_size.
1302  k_min = 0
1303  DO k = k_use - 1, 1, -1
1304  IF (this%delta_a_len(k) .le. step_size) THEN
1305  k_min = k
1306  EXIT
1307  END IF
1308  END DO
1309 
1310 ! Set up the quadratic
1311  ysq = dot_product(this%delta_a(:,k_min), this%delta_a(:,k_min))
1312  zsq = dot_product(this%delta_a(:,k_min + 1), &
1313  & this%delta_a(:,k_min + 1))
1314  ydz = dot_product(this%delta_a(:,k_min), &
1315  & this%delta_a(:,k_min + 1))
1316 
1317 ! Special case when (k_min .eq. 0) - reconstruction arrays contain SD step, but
1318 ! use the zero step for (kmin .eq. 0)
1319  IF (k_min .eq. 0) THEN
1320  ysq = 0.0
1321  ydz = 0.0
1322  END IF
1323 
1324  aa = ysq + zsq - 2.0*ydz
1325  bb = -2.0*ysq + 2.0*ydz
1326  cc = ysq - step_size**2.0
1327 
1328  discriminant = bb**2.0 - 4.0*aa*cc
1329 
1330  IF (discriminant .le. 0.0) THEN
1331  xx = 0.0
1332  ELSE
1333  xx = (-bb + sqrt(discriminant))/(2.0*aa)
1334  END IF
1335 
1336  IF ((xx .lt. 0.0) .or. (xx .gt. 1.0)) THEN
1337  xx = 0.0
1338  END IF
1339 
1340  delta_a = xx*this%delta_a(:,k_min + 1)
1341  IF (k_min .ne. 0) THEN
1342  delta_a = this%delta_a(:,k_min)*(1.0 - xx) + delta_a
1343  END IF
1344 
1345  reconstruction_seg_step = sqrt(dot_product(delta_a, delta_a))
1346  ELSE
1347  delta_a = this%delta_a(:,k_use)
1348  reconstruction_seg_step = this%delta_a_len(k_use)
1349  END IF
1350 
1351  CALL profiler_set_stop_time('reconstruction_seg_step', start_time)
1352 
1353  END FUNCTION
1354 
1355 !-------------------------------------------------------------------------------
1366 !-------------------------------------------------------------------------------
1367  FUNCTION reconstruction_lm_step(this, k_use, step_size, delta_a)
1368 
1369  IMPLICIT NONE
1370 
1371 ! Declare Arguments
1372  REAL (rprec) :: reconstruction_lm_step
1373  TYPE (reconstruction_class), INTENT(inout) :: this
1374  INTEGER, INTENT(in) :: k_use
1375  REAL (rprec), INTENT(in) :: step_size
1376  REAL (rprec), DIMENSION(:), INTENT(out) :: delta_a
1377 
1378 ! local variables
1379  REAL (rprec) :: lambda
1380  REAL (rprec), DIMENSION(:), ALLOCATABLE :: utdote
1381  REAL (rprec) :: start_time
1382 
1383 ! Start of executable code
1384  start_time = profiler_get_start_time()
1385 
1386  this%lm_ratio = 0.0
1387 
1388  IF ((step_size .ge. 0.0) .and. &
1389  & (this%delta_a_len(k_use) .gt. step_size)) THEN
1390 
1391 ! Define a default value of lambda in case the root find does not work well.
1392  lambda = &
1393  & (this%j_svd_w(1)*this%delta_a_len(k_use)/step_size)**2.0
1394 
1395  ALLOCATE(utdote(SIZE(this%e, 1)))
1396  utdote = matmul(this%e(:,this%current_step), this%j_svd_u)
1397 
1398 ! Find the L-M parameter lambda that corresponds to a step length of step_size.
1399  lambda = reconstruction_lm_rootfind(this, k_use, lambda, &
1400  & step_size, utdote)
1401 
1402 ! Find the step
1403  utdote(1:k_use) = utdote(1:k_use)*this%j_svd_w(1:k_use) &
1404  & / (this%j_svd_w(1:k_use)**2.0 + lambda)
1405 
1406  delta_a = matmul(utdote(1:k_use), this%j_svd_vt(1:k_use,:))
1407  reconstruction_lm_step = sqrt(dot_product(delta_a, delta_a))
1408 
1409  this%lm_ratio(1:k_use,1) = lambda
1410  this%lm_ratio(1:k_use,2) = this%j_svd_w(1:k_use)**2.0
1411  this%lm_ratio(1:k_use,3) = this%lm_ratio(1:k_use,1) &
1412  & / this%lm_ratio(1:k_use,2)
1413 
1414  DEALLOCATE(utdote)
1415  ELSE
1416  delta_a = this%delta_a(:,k_use)
1417  reconstruction_lm_step = this%delta_a_len(k_use)
1418  END IF
1419 
1420  CALL profiler_set_stop_time('reconstruction_lm_step', start_time)
1421 
1422  END FUNCTION
1423 
1424 !-------------------------------------------------------------------------------
1436 !-------------------------------------------------------------------------------
1437  FUNCTION reconstruction_lm_rootfind(this, k_use, lambda_default, &
1438  & step_size, utdote)
1439 
1440  IMPLICIT NONE
1441 
1442 ! Declare Arguments
1443  REAL (rprec) :: reconstruction_lm_rootfind
1444  TYPE (reconstruction_class), INTENT(inout) :: this
1445  INTEGER, INTENT(in) :: k_use
1446  REAL (rprec), INTENT(in) :: lambda_default
1447  REAL (rprec), INTENT(in) :: step_size
1448  REAL (rprec), DIMENSION(:), INTENT(in) :: utdote
1449 
1450 ! local variables
1451  INTEGER :: i
1452  REAL (rprec), DIMENSION(:), ALLOCATABLE :: f_sqrd
1453  REAL (rprec) :: step_size_sqrd
1454  REAL (rprec) :: f_0
1455  REAL (rprec) :: f_2
1456  REAL (rprec) :: lambda_0
1457  REAL (rprec) :: lambda_2
1458  REAL (rprec) :: lambda_try
1459  REAL (rprec) :: lambda_new
1460  REAL (rprec) :: f_try, f_new
1461  REAL (rprec) :: s
1462  REAL (rprec) :: start_time
1463 
1464 ! local parameters
1465  REAL (rprec), PARAMETER :: d_lambda_min = 1.0e-6_rprec
1466 
1467 ! Start of executable code
1468  start_time = profiler_get_start_time()
1469 
1470 ! Set up stuff for the function reconstruction_lm_function for root find.
1471  ALLOCATE(f_sqrd(k_use))
1472 
1473 ! Put (U^t*e)^2 in the function.
1474  f_sqrd = utdote(1:k_use)**2.0
1475  step_size_sqrd = step_size**2.0
1476 
1477  reconstruction_lm_rootfind = lambda_default
1478 
1479 ! Find the bracketing values of lambda. Look for a small value of lambda, to
1480 ! give a positive function value. f_0 should be greater than zero, as the step
1481 ! size at k_use is larger than step_size.
1482  lambda_0 = 0.0
1483  f_0 = reconstruction_lm_function(this, f_sqrd, step_size_sqrd, &
1484  & lambda_0)
1485  IF (f_0 .lt. 0.0) THEN
1486  CALL err_warn('reconstruction_lm_rootfind: f_0 < 0.0')
1487  DEALLOCATE(f_sqrd)
1488  CALL profiler_set_stop_time('reconstruction_lm_rootfind', &
1489  & start_time)
1490  RETURN
1491  END IF
1492 
1493 ! Look for a large value of lambda, to give a negative function value.
1494  lambda_2 = this%j_svd_w(1)**2.0
1495  f_2 = reconstruction_lm_function(this, f_sqrd, step_size_sqrd, &
1496  & lambda_2)
1497  DO i = 1, 20
1498  IF (f_2 .le. 0.0) THEN
1499  EXIT
1500  ELSE
1501  lambda_2 = 4.0*lambda_2
1502  f_2 = reconstruction_lm_function(this, f_sqrd, &
1503  & step_size_sqrd, &
1504  & lambda_2)
1505  END IF
1506  END DO
1507  IF (f_2 .gt. 0.0) THEN
1508  CALL err_warn('reconstruction_lm_rootfind: f_2 > 0.0')
1509  DEALLOCATE(f_sqrd)
1510  CALL profiler_set_stop_time('reconstruction_lm_rootfind', &
1511  & start_time)
1512  RETURN
1513  END IF
1514 
1515 ! Check that f_0 and f_2 actually bracket.
1516  IF (f_0*f_2 .gt. 0.0) THEN
1517  CALL err_warn('reconstruction_lm_rootfind: fail to bracket')
1518  DEALLOCATE(f_sqrd)
1519 
1520  CALL profiler_set_stop_time('reconstruction_lm_rootfind', &
1521  & start_time)
1522  RETURN
1523  END IF
1524 
1525 ! Find the root using the Ridder method. Numerical Recipes.
1526  reconstruction_lm_rootfind = -9.99e99_rprec
1527  DO i = 1, 99
1528  lambda_try = 0.5*(lambda_0 + lambda_2)
1529  f_try = reconstruction_lm_function(this, f_sqrd, &
1530  & step_size_sqrd, &
1531  & lambda_try)
1532  s = sqrt(f_try**2.0 - f_0*f_2)
1533  IF (s .eq. 0.0) THEN
1534  EXIT
1535  END IF
1536 
1537  lambda_new = lambda_try + (lambda_try - lambda_0)* &
1538  & (sign(1.0_rprec, f_0 - f_2)*f_try/s)
1539 
1540  IF (abs(lambda_new - reconstruction_lm_rootfind) .le. &
1541  & d_lambda_min) THEN
1542  EXIT
1543  END IF
1544 
1545  reconstruction_lm_rootfind = lambda_new
1546  f_new = reconstruction_lm_function(this, f_sqrd, &
1547  & step_size_sqrd, &
1548  & lambda_new)
1549  IF (f_new .eq. 0.0) THEN
1550  EXIT
1551  ELSE IF (sign(f_try, f_new) .ne. f_try) THEN
1552  lambda_0 = lambda_try
1553  f_0 = f_try
1554  lambda_2 = lambda_new
1555  f_2 = f_new
1556  ELSE IF (sign(f_0, f_new) .ne. f_0) THEN
1557  lambda_2 = lambda_new
1558  f_2 = f_new
1559  ELSE IF (sign(f_2, f_new) .ne. f_2) THEN
1560  lambda_0 = lambda_new
1561  f_0 = f_new
1562  ELSE
1563  stop 'reconstruction_lm_rootfind: should never get here'
1564  END IF
1565 
1566  IF (abs(lambda_0 - lambda_2) .le. d_lambda_min) THEN
1567  EXIT
1568  END IF
1569  END DO
1570 
1571  DEALLOCATE(f_sqrd)
1572 
1573  CALL profiler_set_stop_time('reconstruction_lm_rootfind', &
1574  & start_time)
1575 
1576  END FUNCTION
1577 
1578 !-------------------------------------------------------------------------------
1589 !-------------------------------------------------------------------------------
1590  FUNCTION reconstruction_lm_function(this, f_sqrd, step_size_sqrd, &
1591  & lambda)
1592 
1593  IMPLICIT NONE
1594 
1595 ! Declare Arguments
1596  REAL (rprec) :: reconstruction_lm_function
1597  TYPE (reconstruction_class), INTENT(inout) :: this
1598  REAL (rprec), DIMENSION(:), INTENT(in) :: f_sqrd
1599  REAL (rprec), INTENT(in) :: step_size_sqrd
1600  REAL (rprec), INTENT(in) :: lambda
1601 
1602 ! local variables
1603  INTEGER :: i
1604  REAL (rprec) :: start_time
1605 
1606 ! Start of executable code
1607  start_time = profiler_get_start_time()
1608 
1609  reconstruction_lm_function = 0.0
1610  DO i = 1, SIZE(f_sqrd)
1611  reconstruction_lm_function = reconstruction_lm_function &
1612  & + f_sqrd(i)*(this%j_svd_w(i)/(this%j_svd_w(i)**2.0 + &
1613  & lambda))**2.0
1614  END DO
1615 
1616  reconstruction_lm_function = reconstruction_lm_function &
1617  & - step_size_sqrd
1618 
1619  CALL profiler_set_stop_time('reconstruction_lm_function', &
1620  & start_time)
1621 
1622  END FUNCTION
1623 
1624 !-------------------------------------------------------------------------------
1643 !-------------------------------------------------------------------------------
1644  FUNCTION reconstruction_step(this, signals, derived_params, &
1645  & locks, a_model, gaussp, params, &
1646  & eq_steps, iou, recon_comm, eq_comm)
1647 
1648  IMPLICIT NONE
1649 
1650 ! Declare Arguments
1651  LOGICAL :: reconstruction_step
1652  TYPE (reconstruction_class), INTENT(inout) :: this
1653  TYPE (signal_pointer), DIMENSION(:), INTENT(inout) :: signals
1654  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: &
1655  & derived_params
1656  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: locks
1657  TYPE (model_class), POINTER :: a_model
1658  TYPE (gaussp_class_pointer), DIMENSION(:), INTENT(inout) :: &
1659  & gaussp
1660  TYPE (param_pointer), DIMENSION(:), INTENT(inout) :: params
1661  INTEGER, INTENT(inout) :: eq_steps
1662  INTEGER, INTENT(in) :: iou
1663  INTEGER, INTENT(in) :: recon_comm
1664  INTEGER, INTENT(in) :: eq_comm
1665 
1666 ! local variables
1667  INTEGER :: k_use
1668  INTEGER :: i
1669  INTEGER :: j
1670  REAL (rprec) :: step_use
1671  REAL (rprec) :: start_time
1672  LOGICAL :: temp
1673  REAL (rprec), DIMENSION(4) :: temp_array
1674  INTEGER :: error
1675 
1676 ! Start of executable code
1677  start_time = profiler_get_start_time()
1678  reconstruction_step = .false.
1679 
1680 ! Evaluate the quasi-Newton step.
1681  k_use = reconstruction_eval_step(this, signals, derived_params, &
1682  & locks, a_model, gaussp, params, &
1683  & eq_steps, iou, recon_comm, &
1684  & eq_comm)
1685  this%num_sv(this%current_step + 1) = k_use
1686 
1687 ! Sync the results of the jacobian and singular value decomposition to the
1688 ! child processes.
1689 #if defined(MPI_OPT)
1690  CALL mpi_bcast(mpi_sync_task, 1, mpi_integer, 0, recon_comm, &
1691  & error)
1692  CALL reconstruction_sync_svd(this, recon_comm)
1693  DO i = 1, SIZE(params)
1694  CALL param_sync_delta(params(i)%p, recon_comm)
1695  END DO
1696 #endif
1697 
1698 ! Set the initial step size.
1699  step_use = this%step_max
1700 
1701  DO i = 1, reconstruction_max_step_try
1702 
1703  IF (reconstruction_try_step(this, signals, derived_params, &
1704  & locks, a_model, gaussp, params, &
1705  & eq_steps, step_use, iou, &
1706  & recon_comm, eq_comm)) THEN
1707 
1708 ! Equilibrium converged
1709 
1710  IF (reconstruction_get_dg2(this) .ge. 0.0) THEN
1711  CALL reconstruction_eval_f(this, derived_params, a_model)
1712 
1713 ! The signal cache holds the next set of good signal values.
1714  DO j = 1, SIZE(signals)
1715  this%last_values(:,j) = signals(j)%p%modeled
1716  END DO
1717 
1718  reconstruction_step = .true.
1719  WRITE (*,*)
1720  WRITE (*,*) ' *** Step succeeded'
1721  WRITE (*,*) ' Result of step'
1722  WRITE (*,1000)
1723  WRITE (*,1001) reconstruction_get_lastg2(this), &
1724  & reconstruction_get_g2(this), &
1725  & reconstruction_get_dg2(this)
1726 
1727  WRITE (iou,*)
1728  WRITE (iou,*) ' *** Step succeeded'
1729  WRITE (iou,*) ' Result of step'
1730  WRITE (iou,1000)
1731  WRITE (iou,1001) reconstruction_get_lastg2(this), &
1732  & reconstruction_get_g2(this), &
1733  & reconstruction_get_dg2(this)
1734  EXIT
1735  ELSE
1736 ! Equilibrium converged but g^2 increased
1737  WRITE (*,*) 'Equilibrium converged but g2 increased'
1738  WRITE (*,1000)
1739  WRITE (*,1001) reconstruction_get_lastg2(this), &
1740  & reconstruction_get_g2(this), &
1741  & reconstruction_get_dg2(this)
1742 ! WRITE (*,1002) MAXLOC(ABS(delta_a))
1743  WRITE (iou,*) 'Equilibrium converged but g2 increased'
1744  WRITE (iou,1000)
1745  WRITE (iou,1001) reconstruction_get_lastg2(this), &
1746  & reconstruction_get_g2(this), &
1747  & reconstruction_get_dg2(this)
1748 ! WRITE (iou,1002) MAXLOC(ABS(delta_a))
1749  END IF
1750 
1751  ELSE
1752 ! Equilibrium did not converge
1753  WRITE (*,*) 'Equilibrium did not converge'
1754  WRITE (iou,*) 'Equilibrium did not converge'
1755  END IF
1756 
1757 ! If this section was reached, the reconstruction step failed. Reset the
1758 ! current step to the last good step.
1759  this%current_step = this%current_step - 1
1760 
1761  IF (i .eq. reconstruction_max_step_try) THEN
1762 ! The signals contain cached values. If the step fails recompute the modeled
1763 ! signals so that values written to the recout match the result file.
1764  temp = reconstruction_eval_e(this, signals, a_model, &
1765  & gaussp, eq_steps, iou, eq_comm)
1766 
1767  WRITE (*,*) 'Reconstruction step failed.'
1768  WRITE (iou,*) 'Reconstruction step failed.'
1769  EXIT
1770  END IF
1771 
1772 ! Retry step with a smaller stepsize.
1773  step_use = step_use/2.0
1774  WRITE (*,*) 'Changing max step size to ', step_use
1775  WRITE (iou,*) 'Changing max step size to ', step_use
1776  END DO
1777 
1778  CALL profiler_set_stop_time('reconstruction_step', start_time)
1779 
1780 1000 FORMAT('g^2 old',7x,'g^2 new',7x,'g^2 diff')
1781 1001 FORMAT(es12.5,2x,es12.5,2x,es12.5)
1782 1002 FORMAT('Change vrnc value for rp number ',i4)
1783 
1784  END FUNCTION
1785 
1786 !-------------------------------------------------------------------------------
1805 !-------------------------------------------------------------------------------
1806  FUNCTION reconstruction_try_step(this, signals, derived_params, &
1807  & locks, a_model, gaussp, params, &
1808  & eq_steps, max_step, iou, &
1809  & recon_comm, eq_comm)
1810 
1811  IMPLICIT NONE
1812 
1813 ! Declare Arguments
1814  LOGICAL :: reconstruction_try_step
1815  TYPE (reconstruction_class), INTENT(inout) :: this
1816  TYPE (signal_pointer), DIMENSION(:), INTENT(inout) :: signals
1817  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: &
1818  & derived_params
1819  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: locks
1820  TYPE (model_class), POINTER :: a_model
1821  TYPE (gaussp_class_pointer), DIMENSION(:), INTENT(inout) :: &
1822  & gaussp
1823  TYPE (param_pointer), DIMENSION(:), INTENT(inout) :: params
1824  INTEGER, INTENT(inout) :: eq_steps
1825  REAL (rprec), INTENT(inout) :: max_step
1826  INTEGER, INTENT(in) :: iou
1827  INTEGER, INTENT(in) :: recon_comm
1828  INTEGER, INTENT(in) :: eq_comm
1829 
1830 ! local variables
1831  INTEGER :: i
1832  INTEGER :: mpi_rank
1833  INTEGER :: mpi_size
1834  INTEGER :: error
1835  LOGICAL :: converged
1836  INTEGER :: k_use
1837  REAL (rprec) :: step_use
1838  REAL (rprec) :: delta_norm
1839  REAL (rprec), DIMENSION(:), ALLOCATABLE :: delta_a
1840  REAL (rprec) :: best_g2
1841  REAL (rprec) :: temp_g2
1842  INTEGER :: best_index
1843  REAL (rprec), DIMENSION(:), ALLOCATABLE :: param_value
1844  REAL (rprec) :: start_time
1845 
1846 ! Start of executable code
1847  start_time = profiler_get_start_time()
1848 
1849  IF (this%use_central) THEN
1850  delta_norm = 2.0
1851  ELSE
1852  delta_norm = 1.0
1853  END IF
1854 
1855 #if defined(MPI_OPT)
1856  CALL mpi_comm_rank(recon_comm, mpi_rank, error)
1857  CALL mpi_comm_size(recon_comm, mpi_size, error)
1858 #else
1859  mpi_rank = 0
1860  mpi_size = 1
1861 #endif
1862 
1863  k_use = this%num_sv(this%current_step + 1)
1864 
1865 ! Change the step to use based on the mpi_rank. This heuristic sub divids the
1866 ! max step size to one step just under half the size. For instance if the max
1867 ! step size is 100 and their are two processes, each process will use the
1868 ! following step sizes.
1869 !
1870 ! Process Step Size
1871 ! 0 100.0
1872 ! 1 75.0
1873 !
1874 ! If the are 4 processes, the step sizes will become.
1875 !
1876 ! Process Step Size
1877 ! 0 100.0
1878 ! 1 87.5
1879 ! 2 75.0
1880 ! 3 62.5
1881 !
1882 ! If none of the processes manages to lower g^2, the max step size is halved in
1883 ! reconstruction and the process repeats on the trucated step size. Step size
1884 ! cannot be increased beyond the max step size due to the check in the step
1885 ! types. See the reconstruction_**_step functions.
1886  step_use = max_step - mpi_rank*max_step/(2*mpi_size)
1887 
1888  ALLOCATE(param_value(SIZE(params)))
1889  ALLOCATE(delta_a(SIZE(params)))
1890 
1891 ! Save the parameter values.
1892  DO i = 1, SIZE(params)
1893  param_value(i) = param_get_value(params(i)%p, a_model)
1894  END DO
1895 
1896 ! Step the parameters using the specified step method.
1897  SELECT CASE (this%step_type)
1898 
1899  CASE (reconstruction_sl_step_type)
1900  step_use = reconstruction_sl_step(this, k_use, step_use, &
1901  & delta_a)
1902  this%exp_g2(this%current_step + 1) = &
1903  & reconstruction_get_exp_g2(this, delta_a)
1904  this%step_size(this%current_step + 1) = step_use
1905  CALL reconstruction_write_step(this, 'sl', iou)
1906 
1907  CASE (reconstruction_lm_step_type)
1908  step_use = reconstruction_lm_step(this, k_use, step_use, &
1909  & delta_a)
1910  this%exp_g2(this%current_step + 1) = &
1911  & reconstruction_get_exp_g2(this, delta_a)
1912  this%step_size(this%current_step + 1) = step_use
1913  CALL reconstruction_write_step(this, 'lm', iou)
1914 
1915  CASE (reconstruction_seg_step_type)
1916  step_use = reconstruction_seg_step(this, k_use, step_use, &
1917  & delta_a)
1918  this%exp_g2(this%current_step + 1) = &
1919  & reconstruction_get_exp_g2(this, delta_a)
1920  this%step_size(this%current_step + 1) = step_use
1921  CALL reconstruction_write_step(this, 'seg', iou)
1922 
1923  END SELECT
1924 
1925 ! Rank zero will attempt to take the largest step. Once the largest step is
1926 ! determined, broadcast to children. This avoids the children taking redundant
1927 ! steps.
1928  IF (mpi_rank .eq. 0) THEN
1929 #if defined(MPI_OPT)
1930  CALL mpi_bcast(mpi_step_task, 1, mpi_integer, 0, recon_comm, &
1931  & error)
1932  CALL mpi_bcast(step_use, 1, mpi_real8, 0, recon_comm, error)
1933  CALL mpi_bcast(eq_steps, 1, mpi_integer, 0, recon_comm, error)
1934 #endif
1935  max_step = step_use
1936  END IF
1937 
1938 ! Set the new parameter values. ABS(params(j)%p%delta) is the normalization
1939 ! factor. This is not vrnc because the increment function can take a non vrnc
1940 ! step size.
1941  WRITE (iou,1000)
1942  DO i = 1, SIZE(params)
1943  CALL param_set_value(params(i)%p, a_model, param_value(i) + &
1944  & delta_a(i)*delta_norm*abs(params(i)%p%recon%delta), &
1945  & eq_comm, this%use_central)
1946  WRITE (iou,1001) delta_a(i), param_get_name(params(i)%p, &
1947  & a_model)
1948  END DO
1949  DO i = 1, SIZE(locks)
1950  CALL param_set_lock_value(locks(i)%p, a_model, eq_comm)
1951  END DO
1952 
1953 ! Reconverge the equilibrium. Increment the current step so
1954 ! reconstruction_eval_e writes the e array to the correct index.
1955  this%current_step = this%current_step + 1
1956  converged = reconstruction_eval_e(this, signals, a_model, &
1957  & gaussp, eq_steps, iou, eq_comm)
1958 
1959 ! Search all the processes for the lowest value of g^2.
1960  best_index = -1
1961 #if defined(MPI_OPT)
1962  IF (mpi_rank .gt. 0) THEN
1963  CALL mpi_ssend(converged, 1, mpi_logical, 0, mpi_rank, &
1964  & recon_comm, error)
1965 
1966 ! If the equilibrium converged send the value of g^2.
1967  IF (converged) THEN
1968  best_g2 = reconstruction_get_g2(this)
1969  WRITE (*,1002) mpi_rank, best_g2, step_use
1970  CALL mpi_ssend(best_g2, 1, mpi_real8, 0, mpi_rank, &
1971  & recon_comm, error)
1972  END IF
1973 
1974  ELSE
1975 #endif
1976  IF (converged) THEN
1977  best_index = 0
1978  best_g2 = reconstruction_get_g2(this)
1979  WRITE (*,1002) mpi_rank, best_g2, step_use
1980  END IF
1981 
1982 #if defined(MPI_OPT)
1983  DO i = 1, mpi_size - 1
1984  CALL mpi_recv(converged, 1, mpi_logical, i, i, recon_comm, &
1985  & mpi_status_ignore, error)
1986 
1987 ! If the equilibrium on that processor converged, receive the value of g^2 and
1988 ! reset the best index if g^2 is lower.
1989  IF (converged) THEN
1990  CALL mpi_recv(temp_g2, 1, mpi_real8, i, i, recon_comm, &
1991  & mpi_status_ignore, error)
1992  IF (temp_g2 .lt. best_g2) THEN
1993  best_index = i
1994  best_g2 = temp_g2
1995  END IF
1996  END IF
1997  END DO
1998  END IF
1999 
2000 ! Now that we identified the best index, broadcast that values to all the
2001 ! processes.
2002  CALL mpi_bcast(best_index, 1, mpi_integer, 0, recon_comm, error)
2003 
2004  IF (mpi_rank .eq. 0) THEN
2005  WRITE (*,1003) best_index
2006  IF (best_index .gt. 0) THEN
2007 ! Recieve the best parameters from the best process and sync the model. Store
2008 ! these into param_value. For the zeroth rank, the old parameters were saved in
2009 ! reconstruction_step.
2010  DO i = 1, SIZE(params)
2011  CALL param_sync_child(params(i)%p, a_model, best_index, &
2012  & recon_comm, eq_comm, &
2013  & this%use_central)
2014  END DO
2015  CALL mpi_recv(this%e(:,this%current_step), SIZE(this%e, 1), &
2016  & mpi_real8, best_index, best_index, recon_comm, &
2017  & mpi_status_ignore, error)
2018  CALL mpi_recv(this%exp_g2(this%current_step), 1, mpi_real8, &
2019  & best_index, best_index, recon_comm, &
2020  & mpi_status_ignore, error)
2021 
2022 ! For some reason, this needs to be before the model is synced or it can dead
2023 ! lock.
2024  DO i = 1, SIZE(signals)
2025  CALL signals(i)%p%sync_child(best_index, recon_comm)
2026  END DO
2027 
2028 ! FIXME: Not sure if the model needs to be updated before the lock values are
2029 ! are set.
2030  DO i = 1, SIZE(locks)
2031  CALL param_set_lock_value(locks(i)%p, a_model, eq_comm)
2032  END DO
2033 
2034 ! Sync the model from the best index.
2035  CALL model_sync_child(a_model, best_index, recon_comm)
2036  a_model%state_flags = model_state_all_off
2037  END IF
2038 
2039 ! Ensure that parent had a chance sync the child equilibrium before it gets
2040 ! reset. See below.
2041  CALL mpi_barrier(recon_comm, error)
2042 
2043 ! If the equilibrium failed to converge reset the parameters and equilibrium.
2044  IF ((best_index .lt. 0) .or. &
2045  & (reconstruction_get_dg2(this) .lt. 0.0)) THEN
2046 ! The step failed, reset the param values so that they remain in sync with the
2047 ! step.
2048  DO i = 1, SIZE(params)
2049  CALL param_set_value(params(i)%p, a_model, &
2050  & param_value(i), eq_comm, &
2051  & this%use_central)
2052  END DO
2053  DO i = 1, SIZE(locks)
2054  CALL param_set_lock_value(locks(i)%p, a_model, eq_comm)
2055  END DO
2056 
2057 ! Reset the model state back to the last step. The state of the last step was
2058 ! saved before the jacobian was computed. This needs to be called after the
2059 ! parameters are reset so that the model is flagged to the converged state.
2060 ! Otherwise the model will try to converge again producing an incorrect wout
2061 ! file.
2062  CALL model_reset_state(a_model)
2063  END IF
2064  ELSE
2065 
2066  IF (mpi_rank .eq. best_index) THEN
2067 ! Send the best result back to the parent process.
2068  DO i = 1, SIZE(params)
2069  CALL param_sync_child(params(i)%p, a_model, best_index, &
2070  & recon_comm, eq_comm, &
2071  & this%use_central)
2072  END DO
2073  CALL mpi_ssend(this%e(:,this%current_step), SIZE(this%e, 1), &
2074  & mpi_real8, 0, mpi_rank, recon_comm, error)
2075  CALL mpi_ssend(this%exp_g2(this%current_step), 1, mpi_real8, &
2076  & 0, mpi_rank, recon_comm, error)
2077 
2078 ! For some reason, this needs to be before the model is synced or it can dead
2079 ! lock.
2080  DO i = 1, SIZE(signals)
2081  CALL signals(i)%p%sync_child(best_index, recon_comm)
2082  END DO
2083 
2084 ! Sync the model from the best index. The barrier ensures that child's wout
2085 ! file was copied before the child resets it. See model sync above.
2086  CALL model_sync_child(a_model, mpi_rank, recon_comm)
2087  END IF
2088 
2089 ! Reset everything in case this step failed. When the next jacobian is
2090 ! computed, this will be resynced.
2091  DO i = 1, SIZE(params)
2092  CALL param_set_value(params(i)%p, a_model, param_value(i), &
2093  & eq_comm, this%use_central)
2094  END DO
2095  DO i = 1, SIZE(locks)
2096  CALL param_set_lock_value(locks(i)%p, a_model, eq_comm)
2097  END DO
2098 
2099 ! Ensure that parent had a chance to sync the child equilibrium before it gets
2100 ! reset. See above.
2101  CALL mpi_barrier(recon_comm, error)
2102  CALL model_reset_state(a_model)
2103  this%current_step = this%current_step - 1
2104 
2105  END IF
2106 #endif
2107 
2108  reconstruction_try_step = best_index .ge. 0
2109 
2110  DEALLOCATE(param_value)
2111  DEALLOCATE(delta_a)
2112 
2113  CALL profiler_set_stop_time('reconstruction_try_step', start_time)
2114 
2115 1000 FORMAT('Normalized Step',2x,'param_name')
2116 1001 FORMAT(es12.5,5x,a)
2117 1002 FORMAT(' Trying rank: ',i3,' g^2: ',es12.5,' step size: ',es12.5)
2118 1003 FORMAT(' Best index ',i3)
2119 
2120  END FUNCTION
2121 
2122 !-------------------------------------------------------------------------------
2134 !-------------------------------------------------------------------------------
2135  SUBROUTINE reconstruction_eval_sem(this, params, signals, &
2136  & derived_params)
2137 
2138  IMPLICIT NONE
2139 
2140 ! Declare Arguments
2141  TYPE (reconstruction_class), INTENT(inout) :: this
2142  TYPE (param_pointer), DIMENSION(:), INTENT(inout) :: params
2143  TYPE (signal_pointer), DIMENSION(:), INTENT(inout) :: signals
2144  TYPE (param_pointer), DIMENSION(:), INTENT(inout) :: &
2145  & derived_params
2146 
2147 ! local variables
2148  INTEGER :: j
2149  INTEGER :: i
2150  REAL (rprec), DIMENSION(:,:), POINTER :: cp, cd
2151  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: j_dot_cp
2152  REAL (rprec) :: delta_norm
2153  REAL (rprec) :: start_time
2154 
2155 ! local parameters
2156  CHARACTER (len=*), PARAMETER :: &
2157  & sub_name = 'reconstruction_eval_sem'
2158 
2159 ! Start of executable code
2160  start_time = profiler_get_start_time()
2161 
2162  IF (this%use_central) THEN
2163  delta_norm = 2.0
2164  ELSE
2165  delta_norm = 1.0
2166  END IF
2167 
2168 ! Since the jacobian that will be used is normalized, the signal covaiance
2169 ! must be normalized as well. This results in an idenity matrix of
2170 ! SIZE(signals) X SIZE(signals).
2171 !
2172 ! C_p^-1 = J^T * C_s^-1 * J -> A^T * I * A = A^T * A (1)
2173 !
2174 ! Equation 3 in Hanson et. al. doi:10.1088/0029-5515/49/7/075031
2175 ! Note: This is the same as the hessian that was computed for the
2176 ! reconstruction. Use that value to avoid a matrix multiplication.
2177  ALLOCATE(cp(SIZE(params),SIZE(params)))
2178  cp = this%hessian
2179 
2180 ! Invert C_p^-1 to get the parameter covariance matrix.
2181  CALL reconstruction_invert_matrix(this, cp, sub_name, 'C_p^-1')
2182 
2183 ! From the parameter covariance matrix, map the reconstructed uncertainty to
2184 ! the derived parameters by the derived jacobian.
2185 !
2186 ! C_d = J' * C_p * J'^T (2)
2187  IF (SIZE(derived_params) .gt. 0) THEN
2188  ALLOCATE(cd(SIZE(derived_params), SIZE(derived_params)))
2189  cd = matmul(this%derived_jacobian, &
2190  & matmul(cp, transpose(this%derived_jacobian)))
2191 
2192  DO j = 1, SIZE(derived_params)
2193  derived_params(j)%p%sigma = sqrt(cd(j,j))
2194  derived_params(j)%p%correlation = cd(j,:)
2195  END DO
2196 
2197  DEALLOCATE(cd)
2198  END IF
2199 
2200 ! Normalize derived parameter correlation matrix.
2201  CALL reconstruction_normalize_correlations(derived_params)
2202 
2203 ! Store the parameter covariance. Since the normalized signal covariance matrix
2204 ! was used, the parameter covariance needs to be denormalized.
2205  DO j = 1, SIZE(params)
2206  params(j)%p%sigma = &
2207  & sqrt(cp(j,j)*(delta_norm*params(j)%p%recon%delta)**2.0)
2208  DO i = 1, SIZE(params)
2209  params(j)%p%correlation(i) = &
2210  & cp(j,i)*delta_norm*params(j)%p%recon%delta * &
2211  & delta_norm*params(i)%p%recon%delta
2212  END DO
2213  END DO
2214 
2215 ! Normalize reconstruction parameter correlation matrix.
2216  CALL reconstruction_normalize_correlations(params)
2217 
2218 ! Compute signal effectiveness.
2219  ALLOCATE(j_dot_cp(SIZE(signals),SIZE(params)))
2220  j_dot_cp = matmul(this%jacobian, cp)
2221 
2222 ! The signal effectiveness for each parameter is stored in the parameter
2223 ! itself. Equation 6 in Hanson et. al. doi:10.1002/ctpp.200900059
2224 !
2225 ! (J * cp)^2 / (sigma_p^2 * sigma_s^2) (3)
2226 !
2227 ! Due to use of the normalized jacobian and normalized parameter covariance
2228 ! matrix, the 1/sigma_s^2 term is found in the J^2 term. An extra Pi term must
2229 ! be multiplied out. In terms of the normalized quantities, becomes
2230 !
2231 ! (A * cp_n)^2 * pi^2 / sigma_p^2 (4)
2232  DO j = 1, SIZE(params)
2233  params(j)%p%recon%sem = (j_dot_cp(:,j)**2.0) * &
2234  & (delta_norm*params(j)%p%recon%delta)**2.0 / &
2235  & params(j)%p%sigma**2.0
2236  END DO
2237 
2238  DEALLOCATE(j_dot_cp, cp)
2239 
2240  CALL profiler_set_stop_time(sub_name, start_time)
2241 
2242  END SUBROUTINE
2243 
2244 !-------------------------------------------------------------------------------
2254 !-------------------------------------------------------------------------------
2255  SUBROUTINE reconstruction_invert_matrix(this, matrix, sub, name)
2256 
2257  IMPLICIT NONE
2258 
2259 ! Declare Arguments
2260  TYPE (reconstruction_class), INTENT(inout) :: this
2261  REAL (rprec), DIMENSION(:,:), POINTER :: matrix
2262  CHARACTER (len=*), INTENT(in) :: sub
2263  CHARACTER (len=*), INTENT(in) :: name
2264 
2265 ! local variables
2266  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: svd_vt
2267  REAL (rprec), DIMENSION(:,:), ALLOCATABLE :: svd_u
2268  REAL (rprec), DIMENSION(:), ALLOCATABLE :: svd_work
2269  INTEGER :: num_svd
2270  INTEGER :: work_size
2271  INTEGER :: i, status
2272  INTEGER :: m, n
2273  REAL (rprec) :: start_time
2274 
2275 ! Start of executable code
2276  start_time = profiler_get_start_time()
2277 
2278  m = SIZE(matrix, 1)
2279  n = SIZE(matrix, 2)
2280 
2281  ALLOCATE(svd_vt(n,n))
2282  ALLOCATE(svd_u(m,m))
2283  ALLOCATE(svd_work(1))
2284 
2285  svd_vt = 0.0
2286  svd_u = 0.0
2287  this%svd_w = 0.0
2288  svd_work = 0.0
2289 
2290 ! Find the optimal work size.
2291  CALL dgesvd('All', 'All', m, n, matrix, m, this%svd_w, svd_u, m, &
2292  & svd_vt, n, svd_work, -1, status)
2293  work_size = int(svd_work(1))
2294  DEALLOCATE(svd_work)
2295  ALLOCATE(svd_work(work_size))
2296 
2297 ! Factor the matrix to M = U * W * V^T
2298  CALL dgesvd('All', 'All', m, n, matrix, m, this%svd_w, svd_u, m, &
2299  & svd_vt, n, svd_work, work_size, status)
2300  CALL assert_eq(0, status, sub // ': dgesvd problem when ' // &
2301  & 'inverting ' // name)
2302 
2303 ! If the matrix is rectangular the inverse matrix will be transposed.
2304  DEALLOCATE(matrix)
2305  ALLOCATE(matrix(n,m))
2306 
2307 ! Invert the nonzero elements of the W matrix.
2308  matrix = 0
2309  DO i = 1, SIZE(this%svd_w)
2310  IF (this%svd_w(i) .gt. 0.0 .and. &
2311  & this%svd_w(i)/this%svd_w(1) .ge. this%cut_inv_svd) THEN
2312  matrix(i,i) = 1.0/this%svd_w(i)
2313  END IF
2314  END DO
2315 
2316 ! Invert the matric M^-1 = V * W^-1 * U^T
2317  matrix = matmul(transpose(svd_vt), matmul(matrix, &
2318  & transpose(svd_u)))
2319  DEALLOCATE(svd_vt, svd_u, svd_work)
2320 
2321  CALL profiler_set_stop_time('reconstruction_invert_matrix', &
2322  & start_time)
2323 
2324  END SUBROUTINE
2325 
2326 !-------------------------------------------------------------------------------
2333 !-------------------------------------------------------------------------------
2334  SUBROUTINE reconstruction_normalize_correlations(params)
2335 
2336  IMPLICIT NONE
2337 
2338 ! Declare Arguments
2339  TYPE (param_pointer), DIMENSION(:), INTENT(inout) :: params
2340 
2341 ! local variables
2342  INTEGER :: j, i
2343  REAL (rprec) :: start_time
2344 
2345 ! Start of executable code
2346  start_time = profiler_get_start_time()
2347 
2348  DO j = 1, SIZE(params)
2349  params(j)%p%correlation = params(j)%p%correlation &
2350  & / params(j)%p%sigma
2351  DO i = 1, SIZE(params)
2352  params(j)%p%correlation(i) &
2353  & = params(j)%p%correlation(i)/params(i)%p%sigma
2354  END DO
2355  END DO
2356 
2357  CALL profiler_set_stop_time( &
2358  & 'reconstruction_normalize_correlations', start_time)
2359 
2360  END SUBROUTINE
2361 
2362 !-------------------------------------------------------------------------------
2371 !-------------------------------------------------------------------------------
2372  SUBROUTINE reconstruction_write(this, iou)
2373 
2374  IMPLICIT NONE
2375 
2376 ! Declare Arguments
2377  TYPE (reconstruction_class), POINTER :: this
2378  INTEGER, INTENT(in) :: iou
2379 
2380 ! local variables
2381  INTEGER :: i
2382  REAL (rprec) :: temp_g2
2383  REAL (rprec) :: start_time
2384 
2385 ! Start of executable code
2386  start_time = profiler_get_start_time()
2387 
2388  WRITE (*,*)
2389  WRITE (iou,*)
2390  WRITE (*,*) ' *** Reconstruction Steps'
2391  WRITE (iou,*) ' *** Reconstruction Steps'
2392  WRITE (*,1000)
2393  WRITE (iou,1000)
2394 
2395  temp_g2 = dot_product(this%e(:,0), this%e(:,0))
2396  WRITE (*,1001) 0, temp_g2
2397  WRITE (iou,1001) 0, temp_g2
2398 
2399  DO i = 1, this%current_step
2400  temp_g2 = dot_product(this%e(:,i), this%e(:,i))
2401  WRITE (*,1001) i, temp_g2, this%exp_g2(i), this%num_sv(i), &
2402  & this%step_size(i)
2403  WRITE (iou,1001) i, temp_g2, this%exp_g2(i), this%num_sv(i), &
2404  & this%step_size(i)
2405  END DO
2406 
2407  CALL profiler_set_stop_time('reconstruction_write', start_time)
2408 
2409 1000 FORMAT('irstep',2x,'g^2',11x,'Expected g^2',4x,'sv ',2x, &
2410  & 'normalized step size')
2411 1001 FORMAT(i6,2(2x,es12.5),2x,i3,2x,es12.5)
2412 
2413  END SUBROUTINE
2414 
2415 !-------------------------------------------------------------------------------
2425 !-------------------------------------------------------------------------------
2426  SUBROUTINE reconstruction_write_step1(this, step_type, iou)
2427 
2428  IMPLICIT NONE
2429 
2430 ! Declare Arguments
2431  TYPE (reconstruction_class), INTENT(in) :: this
2432  CHARACTER (len=*), INTENT(in) :: step_type
2433  INTEGER, INTENT(in) :: iou
2434 
2435 ! local variables
2436  REAL (rprec) :: start_time
2437 
2438 ! Start of executable code
2439  start_time = profiler_get_start_time()
2440 
2441  WRITE (*,*)
2442  WRITE (*,1000) this%num_sv(this%current_step + 1)
2443  WRITE (*,1001) step_type, this%step_size(this%current_step + 1)
2444  WRITE (*,1002) this%exp_g2(this%current_step + 1)
2445 
2446  WRITE (iou,*)
2447  WRITE (iou,1000) this%num_sv(this%current_step + 1)
2448  WRITE (iou,1001) step_type, this%step_size(this%current_step + 1)
2449  WRITE (iou,1002) this%exp_g2(this%current_step + 1)
2450 
2451  CALL profiler_set_stop_time('reconstruction_write_step1', &
2452  & start_time)
2453 
2454 1000 FORMAT('Number of singular values used: ',i3)
2455 1001 FORMAT(a4,': Step size limited to ',es12.5)
2456 1002 FORMAT(' Expected g^2 = ',es12.5)
2457 
2458  END SUBROUTINE
2459 
2460 !-------------------------------------------------------------------------------
2469 !-------------------------------------------------------------------------------
2470  SUBROUTINE reconstruction_write_step2(this, iou)
2471 
2472  IMPLICIT NONE
2473 
2474 ! Declare Arguments
2475  TYPE (reconstruction_class), INTENT(in) :: this
2476  INTEGER, INTENT(in) :: iou
2477 
2478 ! local variables
2479  CHARACTER (len=21) :: jac_format
2480  INTEGER :: i
2481  REAL (rprec) :: start_time
2482 
2483 ! Start of executable code
2484  start_time = profiler_get_start_time()
2485 
2486  WRITE (iou,*)
2487  WRITE (iou,*) ' *** Jacobian'
2488  WRITE (jac_format,1000) SIZE(this%jacobian, 2)
2489  DO i = 1, SIZE(this%jacobian, 1)
2490  WRITE (iou,jac_format) i, this%jacobian(i,:)
2491  END DO
2492 
2493  IF (this%step_type .eq. reconstruction_lm_step_type) THEN
2494  WRITE (iou,*)
2495  WRITE (iou,*) ' *** Levenberg-Marquardt Parameters'
2496  DO i = 1, this%num_sv(this%current_step)
2497  WRITE (iou,1003) i, this%lm_ratio(i,:)
2498  END DO
2499  END IF
2500 
2501  WRITE (iou,*)
2502  WRITE (iou,*) ' *** Jacobian Singular Values'
2503  WRITE (iou,1001)
2504  DO i = 1, SIZE(this%j_svd_w)
2505  WRITE (iou,1002) this%j_svd_w(i), i .le. &
2506  & this%num_sv(this%current_step)
2507  END DO
2508  WRITE (iou,*)
2509 
2510  WRITE (iou,*)
2511  WRITE (iou,*) ' *** Derived Jacobian'
2512  WRITE (jac_format,1000) SIZE(this%derived_jacobian, 2)
2513  DO i = 1, SIZE(this%derived_jacobian, 1)
2514  WRITE (iou,jac_format) i, this%derived_jacobian(i,:)
2515  END DO
2516 
2517  WRITE (iou,*)
2518  WRITE (iou,*) ' *** Parameter Covarance Singular Values'
2519  WRITE (iou,1001)
2520  DO i = 1, SIZE(this%svd_w)
2521  WRITE (iou,1002) this%svd_w(i), &
2522  & this%svd_w(i)/this%svd_w(1) .ge. this%cut_inv_svd
2523  END DO
2524  WRITE (iou,*)
2525 
2526  CALL profiler_set_stop_time('reconstruction_write_step2', &
2527  & start_time)
2528 
2529 1000 FORMAT('(i4,',i3,'(2x,es12.5))')
2530 1001 FORMAT('svd value',5x,'used')
2531 1002 FORMAT(es12.5,2x,l)
2532 1003 FORMAT(i3,' lambda = ',es12.5,' W_i^2 = ',es12.5,' ratio = ', &
2533  & es12.5)
2534 
2535  END SUBROUTINE
2536 
2537 !*******************************************************************************
2538 ! NETCDF SUBROUTINES
2539 !*******************************************************************************
2540 !-------------------------------------------------------------------------------
2551 !-------------------------------------------------------------------------------
2552  SUBROUTINE reconstruction_restart(this, result_ncid, current_step, &
2553  & signals, derived_params, &
2554  & a_model)
2555  USE ezcdf
2556 
2557  IMPLICIT NONE
2558 
2559 ! Declare Arguments
2560  TYPE (reconstruction_class), INTENT(inout) :: this
2561  INTEGER, INTENT(in) :: result_ncid
2562  INTEGER, INTENT(in) :: current_step
2563  TYPE (signal_pointer), DIMENSION(:), INTENT(inout) :: signals
2564  TYPE (param_pointer), DIMENSION(:), INTENT(in) :: derived_params
2565  TYPE (model_class), POINTER :: a_model
2566 
2567 ! Local variables
2568  REAL (rprec), DIMENSION(:), ALLOCATABLE :: temp_model
2569  INTEGER :: i
2570  INTEGER :: temp_var_id
2571  INTEGER :: status
2572  REAL (rprec) :: start_time
2573 
2574 ! Start of executable code
2575  start_time = profiler_get_start_time()
2576 
2577  status = nf_inq_varid(result_ncid, 'signal_model_value', &
2578  & temp_var_id)
2579  CALL assert_eq(status, nf_noerr, nf_strerror(status))
2580  status = nf_get_vara_double(result_ncid, temp_var_id, &
2581  & (/ 1, 1, current_step /), &
2582  & (/ 4, SIZE(signals), 1 /), &
2583  & this%last_values)
2584  CALL assert_eq(status, nf_noerr, nf_strerror(status))
2585 
2586 !$OMP PARALLEL DO
2587 !$OMP& SCHEDULE(DYNAMIC)
2588 !$OMP& DEFAULT(SHARED)
2589 !$OMP& PRIVATE(i)
2590  DO i = 1, this%last_para_signal
2591  this%e(i,this%current_step) = &
2592  & signals(i)%p%get_e(a_model, .false., this%last_values(:,i))
2593  END DO
2594 !$OMP END PARALLEL DO
2595  DO i = this%last_para_signal + 1, SIZE(signals)
2596  this%e(i,this%current_step) = &
2597  & signals(i)%p%get_e(a_model, .false., this%last_values(:,i))
2598  END DO
2599 
2600  CALL reconstruction_eval_f(this, derived_params, a_model)
2601 
2602  CALL profiler_set_stop_time('reconstruction_restart', start_time)
2603 
2604  END SUBROUTINE
2605 
2606 !*******************************************************************************
2607 ! MPI SUBROUTINES
2608 !*******************************************************************************
2609 !-------------------------------------------------------------------------------
2617 !-------------------------------------------------------------------------------
2618  SUBROUTINE reconstruction_sync_state(this, recon_comm)
2619 
2620  IMPLICIT NONE
2621 
2622 ! Declare Arguments
2623  TYPE (reconstruction_class), INTENT(inout) :: this
2624  INTEGER, INTENT(in) :: recon_comm
2625 
2626 #if defined(MPI_OPT)
2627 ! local variables
2628  REAL (rprec) :: start_time
2629  INTEGER :: error
2630 
2631 ! Start of executable code
2632  start_time = profiler_get_start_time()
2633 
2634  CALL mpi_bcast(this%current_step, 1, mpi_integer, 0, recon_comm, &
2635  & error)
2636  CALL mpi_bcast(this%e(:,this%current_step), SIZE(this%e, 1), &
2637  & mpi_real8, 0, recon_comm, error)
2638  CALL mpi_bcast(this%f(:,this%current_step), SIZE(this%f, 1), &
2639  & mpi_real8, 0, recon_comm, error)
2640  CALL mpi_bcast(this%last_values, 4*SIZE(this%e, 1), mpi_real8, 0, &
2641  & recon_comm, error)
2642 
2643  CALL profiler_set_stop_time('reconstruction_sync_state', &
2644  & start_time)
2645 
2646 #endif
2647 
2648  END SUBROUTINE
2649 
2650 !-------------------------------------------------------------------------------
2658 !-------------------------------------------------------------------------------
2659  SUBROUTINE reconstruction_sync_svd(this, recon_comm)
2660 
2661  IMPLICIT NONE
2662 
2663 ! Declare Arguments
2664  TYPE (reconstruction_class), INTENT(inout) :: this
2665  INTEGER, INTENT(in) :: recon_comm
2666 
2667 #if defined(MPI_OPT)
2668 ! local variables
2669  REAL (rprec) :: start_time
2670  INTEGER :: error
2671 
2672 ! Start of executable code
2673  start_time = profiler_get_start_time()
2674 
2675  CALL mpi_bcast(this%j_svd_w, SIZE(this%j_svd_w), mpi_real8, 0, &
2676  & recon_comm, error)
2677  CALL mpi_bcast(this%j_svd_u, SIZE(this%j_svd_u), mpi_real8, 0, &
2678  & recon_comm, error)
2679  CALL mpi_bcast(this%j_svd_vt, SIZE(this%j_svd_vt), mpi_real8, 0, &
2680  & recon_comm, error)
2681 
2682  CALL mpi_bcast(this%delta_a_len, SIZE(this%delta_a_len), &
2683  & mpi_real8, 0, recon_comm, error)
2684  CALL mpi_bcast(this%delta_a, SIZE(this%delta_a), mpi_real8, 0, &
2685  & recon_comm, error)
2686 
2687  CALL mpi_bcast(this%jacobian, SIZE(this%jacobian), mpi_real8, 0, &
2688  & recon_comm, error)
2689  CALL mpi_bcast(this%hessian, SIZE(this%hessian), mpi_real8, 0, &
2690  & recon_comm, error)
2691 
2692 ! Number of signular values has been stored in the next step but the
2693 ! step hasn't been updated yet.
2694  CALL mpi_bcast(this%num_sv(this%current_step + 1), 1, mpi_integer, &
2695  & 0, recon_comm, error)
2696 
2697  CALL profiler_set_stop_time('reconstruction_sync_svd', start_time)
2698 
2699 #endif
2700 
2701  END SUBROUTINE
2702 
2703 !-------------------------------------------------------------------------------
2716 !-------------------------------------------------------------------------------
2717  SUBROUTINE reconstruction_sync_parent(this, params, num_signal, &
2718  & num_derived_params, stride, &
2719  & offset, recon_comm)
2720 
2721  IMPLICIT NONE
2722 
2723 ! Declare Arguments
2724  TYPE (reconstruction_class), INTENT(inout) :: this
2725  TYPE (param_pointer), DIMENSION(:), INTENT(inout) :: params
2726  INTEGER, INTENT(in) :: num_signal
2727  INTEGER, INTENT(in) :: num_derived_params
2728  INTEGER, DIMENSION(:), INTENT(in) :: stride
2729  INTEGER, DIMENSION(:), INTENT(in) :: offset
2730  INTEGER, INTENT(in) :: recon_comm
2731 
2732 #if defined(MPI_OPT)
2733 ! local variables
2734  REAL (rprec) :: start_time
2735  INTEGER :: error
2736  INTEGER :: mpi_rank
2737  INTEGER :: i
2738  INTEGER :: iLow, iHigh
2739 
2740 ! Start of executable code
2741  start_time = profiler_get_start_time()
2742 
2743  CALL mpi_comm_rank(recon_comm, mpi_rank, error)
2744 
2745 ! Need to be careful here with strides and offsets. The mpi uses C style array
2746 ! notation but the data is laid out in memory using fortran array notation.
2747 ! This means that all offsets need to be zero referenced. For the jacobians,
2748 ! the arrays are 2D so all strides and offsets need to be multiplied by the
2749 ! relevant column size.
2750  ilow = offset(mpi_rank + 1) + 1
2751  ihigh = offset(mpi_rank + 1) + stride(mpi_rank + 1)
2752 
2753 ! MPI cannot send and recieve to the same buffer. For the root buffer, use the
2754 ! in place specifer. This makes the recieve buffer act as double duty acts as
2755 ! the send and recieve buffer. Otherwise on the root process, the pointers for
2756 ! the buffers would alias.
2757  IF (mpi_rank .eq. 0) THEN
2758  CALL mpi_gatherv(mpi_in_place, &
2759  & num_signal*stride(mpi_rank + 1), mpi_real8, &
2760  & this%jacobian, num_signal*stride, &
2761  & num_signal*offset, mpi_real8, 0, recon_comm, &
2762  & error)
2763 
2764  CALL mpi_gatherv(mpi_in_place, &
2765  & num_derived_params*stride(mpi_rank + 1), &
2766  & mpi_real8, this%derived_jacobian, &
2767  & num_derived_params*stride, &
2768  & num_derived_params*offset, mpi_real8, 0, &
2769  & recon_comm, error)
2770  ELSE
2771  CALL mpi_gatherv(this%jacobian(:,ilow:ihigh), &
2772  & num_signal*stride(mpi_rank + 1), mpi_real8, &
2773  & this%jacobian, num_signal*stride, &
2774  & num_signal*offset, mpi_real8, 0, recon_comm, &
2775  & error)
2776 
2777  CALL mpi_gatherv(this%derived_jacobian(:,ilow:ihigh), &
2778  & num_derived_params*stride(mpi_rank + 1), &
2779  & mpi_real8, this%derived_jacobian, &
2780  & num_derived_params*stride, &
2781  & num_derived_params*offset, mpi_real8, 0, &
2782  & recon_comm, error)
2783  END IF
2784 
2785 ! Sync the delta used to find the jacobians.
2786  IF (mpi_rank .ne. 0) THEN
2787  DO i = offset(mpi_rank + 1) + 1, offset(mpi_rank + 1) + &
2788  & stride(mpi_rank + 1)
2789  CALL param_send_delta(params(i)%p, i, recon_comm)
2790  END DO
2791  END IF
2792 
2793  IF (mpi_rank .eq. 0) THEN
2794  DO i = offset(1) + stride(1) + 1, SIZE(params)
2795  CALL param_recv_delta(params(i)%p, i, recon_comm)
2796  END DO
2797  END IF
2798 
2799  CALL profiler_set_stop_time('reconstruction_sync_parent', &
2800  & start_time)
2801 
2802 #endif
2803 
2804  END SUBROUTINE
2805 
2806  END MODULE
reconstruction::reconstruction_eval_step
integer function reconstruction_eval_step(this, signals, derived_params, locks, a_model, gaussp, params, eq_steps, iou, recon_comm, eq_comm)
Evaluate the proper quasi-Newton step.
Definition: reconstruction.f:1098
reconstruction::reconstruction_lm_function
real(rprec) function reconstruction_lm_function(this, f_sqrd, step_size_sqrd, lambda)
Levenberg-Marquardt function.
Definition: reconstruction.f:1592
profiler
Defines functions for measuring an tabulating performance of function and subroutine calls....
Definition: profiler.f:13
reconstruction::reconstruction_write
subroutine reconstruction_write(this, iou)
Write out a reconstruction.
Definition: reconstruction.f:2373
reconstruction::reconstruction_sync_svd
subroutine reconstruction_sync_svd(this, recon_comm)
Syncronize svd.
Definition: reconstruction.f:2660
reconstruction::reconstruction_get_exp_dg2
real(rprec) function reconstruction_get_exp_dg2(this, delta_a)
Get the expected dg^2 size.
Definition: reconstruction.f:496
reconstruction::reconstruction_eval_e
logical function reconstruction_eval_e(this, signals, a_model, gaussp, eq_steps, iou, eq_comm)
Evaluate the e vector for a reconstruction step.
Definition: reconstruction.f:670
reconstruction::reconstruction_invert_matrix
subroutine reconstruction_invert_matrix(this, matrix, sub, name)
Invert an M x N matrix.
Definition: reconstruction.f:2256
reconstruction::reconstruction_eval_jacobians
subroutine reconstruction_eval_jacobians(this, signals, derived_params, locks, a_model, gaussp, params, eq_steps, iou, recon_comm, eq_comm)
Evaluate the jacobian.
Definition: reconstruction.f:785
guassian_process::gaussp_class_pointer
Pointer to a gaussian process object. Used for creating arrays of gaussian process pointers....
Definition: guassian_process.f:71
reconstruction::reconstruction_get_exp_g2
real(rprec) function reconstruction_get_exp_g2(this, delta_a)
Get the expected g^2 size.
Definition: reconstruction.f:538
reconstruction::reconstruction_get_k_use
integer function reconstruction_get_k_use(this, num_sv)
Get the number of singular values to use.
Definition: reconstruction.f:399
reconstruction::reconstruction_get_g2
real(rprec) function reconstruction_get_g2(this)
Get the current g^2.
Definition: reconstruction.f:570
reconstruction::reconstruction_get_lastg2
real(rprec) function reconstruction_get_lastg2(this)
Get the last g^2.
Definition: reconstruction.f:600
reconstruction::reconstruction_write_step
Interface for writing out step data.
Definition: reconstruction.f:128
reconstruction
Defines the base class of the type reconstruction_class. This class contains the minimization algorit...
Definition: reconstruction.f:12
reconstruction::reconstruction_restart
subroutine reconstruction_restart(this, result_ncid, current_step, signals, derived_params, a_model)
Restart from result file.
Definition: reconstruction.f:2555
reconstruction::reconstruction_write_step2
subroutine reconstruction_write_step2(this, iou)
Write out a final reconstruction step.
Definition: reconstruction.f:2471
guassian_process
Defines the base class of the type guassian_process_class. The guassian_process contains code to comp...
Definition: guassian_process.f:11
reconstruction::reconstruction_class
Base class containing all the data needed to reconstruct a model.
Definition: reconstruction.f:42
reconstruction::reconstruction_construct
type(reconstruction_class) function, pointer reconstruction_construct(num_steps, num_signals, num_derived_parameters, num_parameters, step_type, step_max, cut_svd, cut_eff, cut_marg_eff, cut_delta_a, cut_dg2, last_para_signal, cut_inv_svd, use_central)
Construct a reconstruction_class object.
Definition: reconstruction.f:174
reconstruction::reconstruction_destruct
subroutine reconstruction_destruct(this)
Deconstruct a reconstruction_class object.
Definition: reconstruction.f:280
reconstruction::reconstruction_seg_step_type
integer, parameter reconstruction_seg_step_type
Segmented path reconstruction step type.
Definition: reconstruction.f:29
profiler::profiler_get_start_time
real(rprec) function profiler_get_start_time()
Gets the start time of profiled function.
Definition: profiler.f:194
reconstruction::reconstruction_sl_step_type
integer, parameter reconstruction_sl_step_type
Straight line reconstruction step type.
Definition: reconstruction.f:25
reconstruction::reconstruction_lm_step
real(rprec) function reconstruction_lm_step(this, k_use, step_size, delta_a)
Take a Levenberg-Marquardt step.
Definition: reconstruction.f:1368
reconstruction::reconstruction_sync_parent
subroutine reconstruction_sync_parent(this, params, num_signal, num_derived_params, stride, offset, recon_comm)
Syncronize child state back to parent.
Definition: reconstruction.f:2720
reconstruction::reconstruction_max_step_try
integer, parameter reconstruction_max_step_try
Maximum number of reconstruction step attemps.
Definition: reconstruction.f:32
reconstruction::reconstruction_eval_sem
subroutine reconstruction_eval_sem(this, params, signals, derived_params)
Evaluate the parameter covariance and signal effectiveness.
Definition: reconstruction.f:2137
reconstruction::reconstruction_sl_step
real(rprec) function reconstruction_sl_step(this, k_use, step_size, delta_a)
Take a straight line step.
Definition: reconstruction.f:1236
reconstruction::reconstruction_eval_f
subroutine reconstruction_eval_f(this, derived_params, a_model)
Evaluate the f vector for a reconstruction step.
Definition: reconstruction.f:738
reconstruction::reconstruction_seg_step
real(rprec) function reconstruction_seg_step(this, k_use, step_size, delta_a)
Take a segmented path step.
Definition: reconstruction.f:1278
reconstruction::reconstruction_lm_rootfind
real(rprec) function reconstruction_lm_rootfind(this, k_use, lambda_default, step_size, utdote)
Root find for Levenberg-Marquardt step.
Definition: reconstruction.f:1439
reconstruction::reconstruction_try_step
logical function reconstruction_try_step(this, signals, derived_params, locks, a_model, gaussp, params, eq_steps, max_step, iou, recon_comm, eq_comm)
Try step size.
Definition: reconstruction.f:1810
reconstruction::reconstruction_get_dg2
real(rprec) function reconstruction_get_dg2(this)
Get the change in g^2.
Definition: reconstruction.f:631
profiler::profiler_set_stop_time
subroutine profiler_set_stop_time(symbol_name, start_time)
Gets the end time of profiled function.
Definition: profiler.f:121
reconstruction::reconstruction_no_step_type
integer, parameter reconstruction_no_step_type
No reconstruction step type.
Definition: reconstruction.f:23
reconstruction::reconstruction_write_step1
subroutine reconstruction_write_step1(this, step_type, iou)
Write out a reconstruction step attempt.
Definition: reconstruction.f:2427
reconstruction::reconstruction_step
logical function reconstruction_step(this, signals, derived_params, locks, a_model, gaussp, params, eq_steps, iou, recon_comm, eq_comm)
Step the reconstruction.
Definition: reconstruction.f:1647
reconstruction::reconstruction_lm_step_type
integer, parameter reconstruction_lm_step_type
Levenberg-Marquardt reconstruction step type.
Definition: reconstruction.f:27
guassian_process::gaussp_set_profile
subroutine gaussp_set_profile(this, a_model)
Set the object and coefficient for an index.
Definition: guassian_process.f:354
reconstruction::reconstruction_normalize_correlations
subroutine reconstruction_normalize_correlations(params)
Normalize a correlation matrix.
Definition: reconstruction.f:2335
reconstruction::reconstruction_sync_state
subroutine reconstruction_sync_state(this, recon_comm)
Syncronize state to children.
Definition: reconstruction.f:2619
Generated on Thu Mar 5 2020 15:49:26 for V3FIT by   doxygen 1.8.17