Cppkpo Procedures

Description

This program calculates the term

which is added to the data vector. In this expression

	is the inverse model covariance matrix
	is the model at the kth iterate
	is the initial model

This procedure is called if the parameter iusepo is set to 1.

Flow

1. The program makes a file connection (33) with the file containing the a priori model, and then reads to model in. The format is the same as that for the regular slowness models read in by Hypit89.

2. The a priori model is subtracted from the existing model (the kth iterate)

3. The resultant vector is multiplied by .

   Notes on the 60 Loop:

   Currently, the option irflag = 2 (full inverse covariance) is not allowed - only straight dampers or layerwise inverse covariance is allowed.

   Case of straight dampers (irflag = 0)

   1. nlnum is set to the number of metablocks in the layer.

   2. For each metablock in the layer:
      1. sigmap and sigmas are set to the dampers for P and S slowness, respectively, for the metablock.

      2. playp and plays are set to the value of the difference vector for the metablock.

      3. The product sigmap*playp is stored in the vector play, and the product sigmas*plays is stored in the vector play (offset by nprm spaces).

   3. The difference vector for all of the variables in the layer now store the products in play.

   Case of layerwise damping (irflag = 1)
   1. The index for the appropriate inverse matrix is determined by comparing the layer number (m) with iclay, which stores the first layer of applicability of a given matrix. Thus, iclay(imat+1) gives the first layer that matrix imat+1 is applicable. As long as m is less than this layer, we retain the same value of imat.

   2. nlnum is set to the number of metablocks in the layer.

   3. For each metablock in the layer:
      1. sigmap and sigmas are set to the dampers for P and S slowness, respectively, for the metablock. NOTE: You shouldn't use this part of the code if idmptp = 2! The damper must be constant on a layerwise basis in order for this procedure to be applicable.

      2. The products of rmm with the appropriate parts of po are performed and stored in playp and plays. Note that rmm is lower triagonal. Thus, if we represent the matrix multiplication as:
         

         The index l of the lower triangle representation of R is computed from

         l = i*(i-1) + k

         when i is greater than k. Because the matrix is symmetric, then the formula

         l = k*(k-1) + i

         will also work, and should be used when i is less than k. When i = k, either can be used.

   4. At this point, the same procedure is done as for irflag = 0.

   At the completion of this routine, the vector po is ready to be added to the data vector (see BUILDG2 procedures for how this is done).

   Authors

   Steve Roecker
   roecks@rpi.edu