az_matrix_util.c

并行解法器,功能强大

  /* do the transpose */

  for (i = 0; i < *n; i++)
    for (j = 0; j < *m; j++) {
      A_index = i + j * *n;
      *(work + w_index++) = *(A + A_index);
    }

  /* copy from "work" back to "A" */

  for (i = 0; i < *m * *n; i++)
    *(A + i) = *(work + i);


  AZ_free((void *) work);

  /* exchange m and n */

  itemp = *m;
  *m    = *n;
  *n    = itemp;

} /* dtrans */

/******************************************************************************/
/******************************************************************************/
/******************************************************************************/

int AZ_get_sym_indx(int iblk, int jblk, int *indx, int *bindx, int *bpntr)

/*******************************************************************************

  Given a block-row and block-column index, return the index to the symmetric
  starting point of the matrix.

  Author:          Scott A. Hutchinson, SNL, 1421
  =======

  Return code:     int, starting index into val of the symmetric portion of the
  ============          matrix.

  Parameter list:
  ===============

  iblk:            Current block row index.

  jblk:            Current block column index.

  indx,
  bindx,
  bpntr:           Arrays used for DMSR and DVBR sparse matrix storage (see
                   file Aztec User's Guide).

*******************************************************************************/

{

  register int i, icount = 0;
  int          itemp, pre_num_nz_blks = 0, num_nz_blks;
  int         *bindx_ptr;

  itemp     = bpntr[jblk];
  bindx_ptr = bindx + itemp;

  /*
   * Count up how many nonzero block precede the current iblk' in the current
   * block row.
   */

  num_nz_blks = bpntr[jblk+1] - itemp;

  for (i = 0; i < num_nz_blks; i++) {
    if (*(bindx_ptr + icount) == iblk) {
      pre_num_nz_blks = icount;
      break;
    }
    else
      icount++;
  }

  return indx[itemp + pre_num_nz_blks];

} /* AZ_get_sym_indx */

/******************************************************************************/
/******************************************************************************/
/******************************************************************************/

extern int AZ_sys_msg_type;

void AZ_print_out(int update_index[], int extern_index[], int update[], 
	int external[], double val[], int indx[], int bindx[], int rpntr[], 
	int cpntr[], int bpntr[], int proc_config[], int choice, int matrix, 
	int N_update, int N_external, int of)
{
/*******************************************************************************

  Print out the matrix in 1 of 3 formats described below.
  starting point of the matrix.

  Author:          Ray Tuminaro, SNL, 9222
  =======

  Return code:     none.
  ============ 

  Parameter list:
  ===============

  update_index,    On input, ordering of update and external locally on this
  extern_index:    processor. For example  'update_index[i]' gives the index
                   location of the block which has the global index 'update[i]'.
                   (AZ_global_mat only).
 
  update:          On input, blks updated on this node (AZ_global_mat only).

  external:        On input, list of external blocks (AZ_global_mat only).

  val,bindx        On input, matrix (MSR or VBR) arrays holding matrix values.
  indx, bnptr,     When using either AZ_input_form or AZ_explicit, these can
  rnptr, cnptr:    be either pre or post-AZ_transform() values depending on what
                   the user wants to see. When using AZ_global_form, these must
                   be post-AZ_transform() values.

  proc_config:     On input, processor information corresponding to:
                      proc_config[AZ_node] = name of this processor
                      proc_config[AZ_N_procs] = # of processors used
 
  choice:          On input, 'choice' determines the output to be printed
		   as described above.
 
  matrix:          On input, type of matrix (AZ_MSR_MATRIX or AZ_VBR_MATRIX).

  N_update:        On input, number of points/blks to be updated on this node.

  N_external:      On input, number of external points/blks needed by this node.

  of:              On input, an offset used with AZ_global_matrix and 
		   AZ_explicit. In particular, a(of,of) is printed for 
                   the matrix element stored as a(0,0).

*******************************************************************************/

   int type, neighbor, cflag;
   int ii,i = 1,j,tj;
   int iblk, jblk, m1, n1, ival, new_iblk, new_jblk;
   MPI_AZRequest request;  /* Message handle */

 
   type            = AZ_sys_msg_type;
   AZ_sys_msg_type =(AZ_sys_msg_type+1-AZ_MSG_TYPE) % AZ_NUM_MSGS +AZ_MSG_TYPE;

   /* Synchronize things so that processor 0 prints first and then */
   /* sends a message to processor 1 so that he prints, etc.      */

   neighbor = proc_config[AZ_node] - 1;
   if ( proc_config[AZ_node] != 0) {
      mdwrap_iread((void *) &i, 0, &neighbor, &type, &request);
      mdwrap_wait((void *) &i, 0, &neighbor, &type, &cflag, &request);
   }
   printf("proc %d:\n",proc_config[AZ_node]);

   if (choice == AZ_input_form ) {
     if ( update != (int *) NULL) {
        printf("  N_update: %5d\n",N_update); printf("  update:");
        AZ_list_print(update, N_update, (double *) NULL , 0);
     }

     if (matrix == AZ_MSR_MATRIX) {
        printf("  bindx: ");
        AZ_list_print(bindx, bindx[N_update], (double *) NULL , 0);

        printf("  val:   ");
        AZ_list_print((int *) NULL , N_update, val , bindx[N_update]);
     }
     else if (matrix == AZ_VBR_MATRIX) {
        printf("  rpntr: ");
        AZ_list_print(rpntr, N_update+1, (double *) NULL , 0);
        if ( cpntr != (int *) NULL ) {
           printf("  cpntr: ");
           AZ_list_print(cpntr, N_update+1+ N_external, (double *) NULL , 0);
        }
        printf("  bpntr: ");
        AZ_list_print(bpntr, N_update+1, (double *) NULL , 0);
        printf("  bindx: ");
        AZ_list_print(bindx, bpntr[N_update], (double *) NULL , 0);
        printf("  indx:  ");
        AZ_list_print(indx, bpntr[N_update]+1, (double *) NULL , 0);
        printf("  val:   ");
        AZ_list_print((int *) NULL, indx[bpntr[N_update]], val, 0);
     }
   }
   else if (choice == AZ_global_mat ) {
     if ( matrix == AZ_MSR_MATRIX) {
        for (i = 0 ; i < N_update; i++ ) {
          ii = update_index[i];
          printf("   a(%d,%d) = %20.13e;\n",update[i]+of,update[i]+of,val[ii]);
          for (j = bindx[ii] ; j < bindx[ii+1] ; j++ ) {
            tj = AZ_find_simple(bindx[j], update_index, N_update, extern_index,
                              N_external,update,external);
            if (tj != -1) 
               printf("   a(%d,%d) = %20.13e;\n",update[i]+of,tj+of,val[j]);
            else (void) fprintf(stderr,"col %d (= bindx[%d]) is undefined\n",
                                tj, j);
          }
        }
     }
     else if (matrix == AZ_VBR_MATRIX) {
        for (iblk= 0; iblk < N_update; iblk++) {
           new_iblk = update_index[iblk];

           m1 = rpntr[new_iblk+1] - rpntr[new_iblk];
 
           /* loop over blocks in the current block-row */
 
           for (ii = bpntr[new_iblk]; ii < bpntr[new_iblk+1]; ii++) {
              new_jblk = AZ_find_simple(bindx[ii], update_index, N_update, 
			 extern_index, N_external,update,external);
              if (new_jblk == -1) {
                 printf("local column %d not found\n",new_jblk);
                 exit(-1);
              }
              jblk = bindx[ii];
              ival =  indx[ii];

              n1 = cpntr[jblk+1] - cpntr[jblk];

              for (i = 0; i < m1; i++) {
                 for (j = 0; j < n1; j++)
                    (void) printf("   a(%d(%d),%d(%d)) = %20.13e;\n",update[iblk]+
				  of,i+of, new_jblk+of, j+of, val[ival+j*m1+i]);
              }
           }
        }
     }
   }
   else if (choice == AZ_explicit) {
     if ( matrix == AZ_MSR_MATRIX) {
        for (i = 0 ; i < N_update; i++ ) {
          if (update == NULL) tj = i+of;
          else tj = update[i] + of;
          printf("   a(%d,%d) = %20.13e;\n",tj,tj,val[i]);
          for (j = bindx[i] ; j < bindx[i+1] ; j++ ) {
               printf("   a(%d,%d) = %20.13e;\n",tj,bindx[j]+of,val[j]);
          }
        }
     }
     else if (matrix == AZ_VBR_MATRIX) {
        for (iblk = 0; iblk < N_update; iblk++) {
           if (update == NULL) tj = iblk+of;
           else tj = update[iblk] + of;

           m1 = rpntr[iblk+1] - rpntr[iblk];
 
           /* loop over blocks in the current block-row */
 
           for (ii = bpntr[iblk]; ii < bpntr[iblk+1]; ii++) {
              jblk = bindx[ii];
              ival =  indx[ii];
              n1 = (indx[ii+1]-ival)/m1;

              for (i = 0; i < m1; i++) {
                 for (j = 0; j < n1; j++)
                    (void) printf("   a(%d(%d),%d(%d)) = %20.13e;\n", tj, 
				  i+of, jblk+of, j+of, val[ival+j*m1+i]);
              }
 
           }
        }
     }
   }
   else (void) fprintf(stderr,"AZ_matrix_out: output choice unknown\n");

   neighbor = proc_config[AZ_node] + 1;
   if ( proc_config[AZ_node] != proc_config[AZ_N_procs] - 1) 
      mdwrap_write((char *) &i, 0, neighbor, type, &cflag);

   i = AZ_gsum_int(i,proc_config);


}
int AZ_find_simple(int k, int *update_index, int N_update, int *extern_index,
   int N_external, int *update, int *external)
{
   int i;

   if (k < N_update) {
      for (i = 0 ; i < N_update ; i++ ) 
         if ( update_index[i] == k) return(update[i]);
   }
   else {
      for (i = 0 ; i < N_external ; i++ ) 
         if ( extern_index[i] == k) return(external[i]);
   }

   return(-1);
}
void AZ_list_print(int ivec[] , int length1, double dvec[],int length2)

{
   int i, count = 0;

   if (ivec != (int *) NULL) {
      for (i = 0 ; i < length1 ; i++ ) {
         printf("%8d ",ivec[i]); count += 8;
         if (count > 50) { count = 0; printf("\n         "); }
      }
   }
   else if (dvec != (double *) NULL) {
      for (i = 0 ; i < length1 ; i++ ) {
         printf("%8.1e ",dvec[i]); count += 8;
         if (count > 50) { count = 0; printf("\n         "); }
      }
      if (length2 != 0) { 
         printf("      -- "); count += 8;
         if (count > 50) { count = 0; printf("\n         "); }
      }

      for (i = length1+1 ; i < length2 ; i++ ) {
         printf("%8.1e ",dvec[i]); count += 8;
         if (count > 50) { count = 0; printf("\n         "); }
      }
   }
   printf("\n");
}

/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
 
int AZ_compute_max_nz_per_row(AZ_MATRIX *Amat, int N, int Nb, int *largest_band)
{
/*******************************************************************************
 
  Compute (and return) the maximum number of nonzeros in 
  any matrix row. Also compute the largest band in the matrix 
  and return its value in '*largest_band'.
 
  Author:          Ray Tuminaro, SNL, 9222
  =======
 
  Return code:     int


  Parameter list:
  ===============
 
  matrix_type:     On input, indicates whether we have an MSR or VBR matrix.

  N:               On input, indicates order of matrix system.

  Nb:              On input, indicates number of blocks for a VBR matrix.

  bindx,bpntr      On input, matrix for which the maximum number of nonzeros
  cpntr:           in a row will be computed.

  largest_band:    On output, largest band in the matrix

*******************************************************************************/

   int i,j,largest = 0;
   int col, col_ptr;
   int kk,left_most,right_most;
   int matrix_type, *bindx, *bpntr, *cpntr;
 
   matrix_type = Amat->matrix_type;
   bindx       = Amat->bindx;
   *largest_band = -1;
   if (matrix_type == AZ_MSR_MATRIX) {
      for (i = 0 ; i < N; i++ ) {
         left_most  = i;
         right_most = i;
         j = bindx[i+1]-bindx[i];
         if (largest < j) largest = j;
         for (kk = bindx[i] ; kk < bindx[i+1] ; kk++ ) {
            if ( bindx[kk] < left_most  ) left_most  = bindx[kk];
            if ( bindx[kk] > right_most ) right_most = bindx[kk];
         }
         if (*largest_band <= right_most - left_most)
            *largest_band = right_most - left_most + 1;
      }
   }
   else if (matrix_type == AZ_VBR_MATRIX) {
      bpntr = Amat->bpntr;
      cpntr = Amat->cpntr;
      for  (i = 0; i < Nb ; i++) {
         j = 0;
         left_most  = cpntr[i];
         right_most = cpntr[i];
         for (col_ptr = bpntr[i]; col_ptr < bpntr[i+1]; col_ptr++ ) {
            col  = bindx[col_ptr];
            if (cpntr[col]   < left_most ) left_most  = cpntr[col];
            if (cpntr[col+1] > right_most) right_most = cpntr[col+1] -1;
            j   += (cpntr[col+1] - cpntr[col]);
         }
         if (*largest_band <= right_most - left_most) 
            *largest_band = right_most - left_most + 1;
         if (j > largest) largest = j;
      }
   }
   else {
      if ( (Amat->largest_band == -1) || (Amat->max_per_row == -1))
          AZ_matfree_Nnzs(Amat);
      *largest_band = Amat->largest_band;
      largest       = Amat->max_per_row - 1;
   }
   (*largest_band)++;
   largest++;
   return(largest);
}

/****************************************************************/
/****************************************************************/
/****************************************************************/
void AZ_check_block_sizes(int bindx[], int cpntr[], int Nrows, 
			  int *new_block)
{
/*
 * This routine is designed for mem_efficient_msr_2_vbr conversion.
 * It takes cpntr[] which is composed of a set of subsequences
 * (a subsequence is set of consecutive cpntr[] values that are 
 * the same) and splits these subsequences if they do not corrspond
 * to blocks in the MSR matrix (bindx,val). 
 *
 * Parameters
 * ==========
 *    bindx,          On input, an msr-LIKE matrix for which we wish
 *                    to discover the block structure. NONZEROS are
 *                    stored consecutively (row-wise) including the 
 *                    diagonal. The last column number in each row
 *                    is encoded as a negative (=  -1 - column).
 *
 *    cpntr[]         On input, cpntr[0] = 0. If the sparsity pattern 
 *                    of row i and row i-1 are identical, cpntr[i] = 
 *                    cpntr[i-1]. Otherwise, cpntr[i] = cpntr[i-1]+1.
 *                    On output, cpntr[0] = 0. If row i and row i-1 can
 *                    be grouped into a block, cpntr[i] = cpntr[i-1].
 *                    Otherwise, cpntr[i] = cpntr[i-1]+1.