az_matrix_util.c

并行解法器,功能强大

 *
 *    Nrows           On input, number of rows in the matrix.
 * 
 *    *new_block      On input, number of unique values in cpntr[].
 *                    On output, number of unique values in the new cpntr[].
 */
 
   int start_blk_ptr, end_blk_ptr, start_blk_column, end_blk_column;
   int blk_name0, blk_name1, blk_name2, blk_name3;
   int prev_col, current_col, start_row, end_row;
   int row, tptr, max_col, last_nz_in_row, N_blks, current;

   max_col = Nrows-1; 

   if (Nrows == 0) return;

   start_blk_ptr = 0;
   row = 0;

   while (row < Nrows) {          /* for each row */
      start_row = start_blk_ptr;
      last_nz_in_row = 0;

      while(!last_nz_in_row) {    /* for each nonzero within a row */
         start_blk_column = bindx[start_blk_ptr];
         if (start_blk_column < 0) { 
            start_blk_column = -1 - start_blk_column; 
            last_nz_in_row = 1;
         }
         prev_col = start_blk_column;
         blk_name1   = cpntr[start_blk_column];

         end_blk_ptr = start_blk_ptr + 1;

         /* Find the end of the current block. If the column numbers */
         /* corresponding to this block are not all found (i.e.      */
         /* current_col != prev_col+1) cut the block off.            */

         while (!last_nz_in_row) {
            current_col = bindx[end_blk_ptr];
            if (current_col < 0) { 
                current_col = -1 - current_col; 
		last_nz_in_row = 1;
            }
            if (current_col != prev_col+1) break;
            if (cpntr[current_col] != blk_name1) break;
            end_blk_ptr++;
            prev_col = current_col;

         }
         end_blk_ptr--;

         /* check to see if the block name associated with the  */
         /* beginning and end of this block are different from  */
         /* the previous and next blocks. If they are the same  */
         /* create a new block name for this block.             */

         end_blk_column = bindx[end_blk_ptr];
         if (end_blk_column < 0) {
            end_blk_column = -1 - end_blk_column;
            last_nz_in_row = 1;
         }
         else last_nz_in_row = 0;

         blk_name2 = cpntr[end_blk_column];

         blk_name0 = -10; blk_name3 = -10;
         if (start_blk_column !=       0) blk_name0 = cpntr[start_blk_column-1];
         if (end_blk_column != max_col) blk_name3 = cpntr[end_blk_column+1];

         if (blk_name1 == blk_name0) {
            for ( tptr = start_blk_column; tptr <= end_blk_column; tptr++) 
               cpntr[tptr] = *new_block;
            (*new_block)++;
         }
         else if (blk_name2 == blk_name3) {
            for ( tptr = start_blk_column; tptr <= end_blk_column; tptr++) 
               cpntr[tptr] = *new_block;
            (*new_block)++;
         }
         start_blk_ptr = end_blk_ptr + 1;
      }
      end_row = end_blk_ptr;

      /* look for the next row which does not have the same */
      /* sparsity pattern as the previous row.              */

      row++;
      while ((row < Nrows) && (cpntr[row-1] == cpntr[row])) {
         row++;
         start_blk_ptr += (end_row - start_row + 1);
      }
   }

   /* Number each of the subsequences in cpntr (a subsequence  */
   /* is series of consecutive cpntr values that are the same) */
   /* uniquely starting from 0 and counting up.                */

   N_blks = 0;
   current = cpntr[0];
   cpntr[0] = N_blks;
   for (row = 1 ; row < Nrows ; row++ ) {
      if (cpntr[row] != current) {
         N_blks++;
         current = cpntr[row];
      }
      cpntr[row] = N_blks;
   }
   *new_block = N_blks;
}
void AZ_capture_matrix(double val[], int indx[], int bindx[], int rpntr[],
                      int cpntr[], int bpntr[], int proc_config[],
                      int data_org[], double b[])

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

  Routine to capture matrix in a form acceptable to MATLAB.  Prints matrix out
  in i,j,a(i,j) format.  Currently implemented only for one processor.  Based
  on the routine AZ_output_matrix.
  Test to see if we should capture matrix, rhs and partitioning info
  in an ASCII data file.  If the file "AZ_write_matrix_now" exists in
  the current working directory, then the files 

  - AZ_capture_matrix.dat
  - AZ_capture_rhs.dat
  - AZ_capture_partition.dat (VBR only)
  
  will be appended with the current matrix in (i,j,val) format, the
  current RHS and the current partition information.  The existence
  of "AZ_write_matrix_now" is check each time.  Thus, capturing can
  be turned on and off at will during the run of a simulation.


  Author:          Michael A. Heroux, 9222, SNL.
  =======

  Return code:     void
  ============

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

  val:             Array containing the nonzero entries of the matrix.

  indx,
  bindx,
  rpntr,
  cpntr,
  bpntr:           Arrays used for DMSR and DVBR sparse matrix storage.

  proc_config:     Machine configuration.  proc_config[AZ_node] is the node
                   number.  proc_config[AZ_N_procs] is the number of processors.

  data_org:        Array containing information on the distribution of the
                   matrix to this processor as well as communication parameters
                   (see file params.txt).
  b:               Right hand side values.

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

{

  /* local variables */

  int  iblk_row, i, j, ib1, ib2, n1, jblk, m1, ipoint, jpoint;
  int  ival = 0;
  int  num_total_nonzeros;
  int  num_total_nodes, N_external_nodes, num_total_equations = 0;
  int  Proc, Num_Proc;

  /********** execution begins **********/
  { FILE *AZ_capture_flag;
  AZ_capture_flag = fopen("AZ_write_matrix_now","r");
  if(AZ_capture_flag)
    {

      Proc               = proc_config[AZ_node];
      Num_Proc           = proc_config[AZ_N_procs];
      if (Num_Proc != 1)
	printf("\nMatrix Capture Routine currently only works for one PE\n");
      else {

	AZ_print_sync_start(Proc, AZ_TRUE, proc_config);

	{ FILE *AZ_capt_matrix, *AZ_capture_rhs, *AZ_capture_partition;

	if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) {
	  num_total_nodes    = data_org[AZ_N_int_blk]+data_org[AZ_N_bord_blk];
	  N_external_nodes = data_org[AZ_N_ext_blk];
	  num_total_equations = rpntr[num_total_nodes];
	  num_total_nonzeros = indx[bpntr[num_total_nodes]];

     /***** Print out the VBR partitioning information for the matrix *****/

	  AZ_capture_partition = fopen("AZ_capture_partition.dat","a");
	  
	  fprintf(AZ_capture_partition, "Start of partition\n");

	  for (i = 0; i < num_total_nodes + 1; i++)
	    fprintf(AZ_capture_partition,"%d\n", rpntr[i]);

	  fclose(AZ_capture_partition);

     /***** Print out the VBR i,j,a(i,j) information for the matrix *****/

	  AZ_capt_matrix = fopen("AZ_capture_matrix.dat","a");
	  fprintf(AZ_capt_matrix, "Start of VBR matrix\n");
	  fprintf(AZ_capt_matrix, "%d %d\n", 
		  num_total_equations, num_total_nonzeros);

	  /* loop over block rows */

	  for (iblk_row = 0; iblk_row < num_total_nodes; iblk_row++) {

	    /* the starting point row index of the current block */

	    ib1 = rpntr[iblk_row];      

	    /* number of rows in the current row block */

	    m1 = rpntr[iblk_row+1] - ib1;

	    /* starting index of current row block */

	    ival = indx[bpntr[iblk_row]];

	    /* loop over all the blocks in the current block-row */

	    for (j = bpntr[iblk_row]; j < bpntr[iblk_row+1]; j++) {
	      jblk = bindx[j];

	      /* the starting point column index of the current block */

	      ib2 = cpntr[jblk];

	      /* number of columns in the current block */

	      n1 = cpntr[jblk+1] - ib2;

	      for (jpoint = 0; jpoint < n1; jpoint++)
		for (ipoint = 0; ipoint < m1; ipoint++) {
		  fprintf(AZ_capt_matrix,"%d %d %22.16e\n", 
			  ib1+ipoint+1, 
			  ib2+jpoint+1, 
			  val[ival+jpoint*m1+ipoint]);
	  
		}

	      ival += m1*n1;
	    }
	  }
	  fclose(AZ_capt_matrix);
	}

	if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {
	  num_total_equations = data_org[AZ_N_internal]+data_org[AZ_N_border];

	  N_external_nodes = data_org[AZ_N_external];
	  num_total_nonzeros = bindx[num_total_equations]-1;

     /***** Print out the MSR i,j,a(i,j) information for the matrix *****/

	  AZ_capt_matrix = fopen("AZ_capture_matrix.dat","a");
	  fprintf(AZ_capt_matrix, "Start of MSR matrix\n");
	  fprintf(AZ_capt_matrix, "%d %d\n", 
		  num_total_equations, num_total_nonzeros);
	  for (i = 0; i < num_total_equations; i++) {
	    fprintf(AZ_capt_matrix,"%d %d %22.16e\n", i+1, i+1, val[i]);
	    for (j = bindx[i]; j < bindx[i+1]; j++ )
	      fprintf(AZ_capt_matrix,"%d %d %22.16e\n", 
		      i+1, bindx[j]+1, val[j]);
	  }
	  fclose(AZ_capt_matrix);
	}

     /***** Print out the RHS information for the matrix *****/

	AZ_capture_rhs = fopen("AZ_capture_rhs.dat","a");
	fprintf(AZ_capture_rhs, "Start of RHS\n");
	for (i = 0; i < num_total_equations; i++)
	  fprintf(AZ_capture_rhs,"%22.16e\n", b[i]);
	fclose(AZ_capture_rhs);

	AZ_print_sync_end(proc_config, AZ_TRUE);
	}
      }
      fclose(AZ_capture_flag);
    }
  }
} /* AZ_capture_matrix */
     

   
struct submat_struct {
   int Nsub_rows, *sub_rows, Nsub_cols, *sub_cols;
};

int AZ_subMSR_getrow(int[], double[], int[], AZ_MATRIX *, int, int[], int);
void AZ_subMSR_matvec_mult (double *, double *, struct AZ_MATRIX_STRUCT *,
	int *);
int AZ_blockMSR_getrow(int[], double[], int[], AZ_MATRIX *, int, int[], 
	int);
void AZ_blockMSR_matvec_mult (double *, double *, struct AZ_MATRIX_STRUCT *,
	 int *);


AZ_MATRIX *AZ_submatrix_create(AZ_MATRIX *Amat, int Nsub_rows, int sub_rows[], 
				int Nsub_cols, int sub_cols[], int proc_config[])
{
	int i;
	
	AZ_MATRIX *sub_matrix;
	struct submat_struct *new_dataptr;

	
	sub_matrix = AZ_matrix_create(Nsub_rows);
	
	new_dataptr = (struct submat_struct *)malloc(sizeof(struct submat_struct));
	
	new_dataptr->Nsub_rows = Nsub_rows;
	new_dataptr->Nsub_cols = Nsub_cols;
	new_dataptr->sub_rows = (int *) malloc(Nsub_rows * sizeof(int));
	new_dataptr->sub_cols = (int *) malloc(Nsub_cols * sizeof(int));

	if ((new_dataptr->sub_rows == NULL) || (new_dataptr->sub_cols == NULL)) {
		printf("couldn't allocate memory for submatrix rows or cols\n");
		exit(-1);
	}

	for (i=0; i < Nsub_rows; i++)
		new_dataptr->sub_rows[i]=sub_rows[i];
	for (i = 0; i < Nsub_cols; i++)
		new_dataptr->sub_cols[i]=sub_cols[i];
	
	sub_matrix->bindx = Amat->bindx;
	sub_matrix->val = Amat->val;

	AZ_set_MATFREE(sub_matrix, new_dataptr, AZ_subMSR_matvec_mult);
	AZ_set_MATFREE_getrow(sub_matrix, new_dataptr, AZ_subMSR_getrow, NULL, 0, proc_config);

	return(sub_matrix);
}

void AZ_submatrix_destroy(AZ_MATRIX **submat)
{
	int *sub_rows, *sub_cols;
	struct submat_struct *data;

	data = (struct submat_struct *) AZ_get_matvec_data(*submat); 

	if (data != NULL) {
		sub_rows = data->sub_rows;
		sub_cols = data->sub_cols;
		
		free(sub_rows);
		free(sub_cols);
		free(data);
	}

	AZ_matrix_destroy(submat);
}
	
int  AZ_subMSR_getrow(int columns[], double values[], int row_lengths[],
											AZ_MATRIX *Amat, int N_requested_rows,
											int requested_rows[], int allocated_space)
{ 
	int    *bindx, i, j, count = 0, fullrow, count_row;
	int     Nsub_rows, Nsub_cols, *sub_rows, *sub_cols, subrow;
	int     fullcol, subcol;
	double *val;
	struct submat_struct *dataptr;
	
	bindx = Amat->bindx;
	val   = Amat->val;
	dataptr = (struct submat_struct *) AZ_get_matvec_data(Amat);
	sub_rows=dataptr->sub_rows;
	sub_cols=dataptr->sub_cols;
	Nsub_rows=dataptr->Nsub_rows;
	Nsub_cols=dataptr->Nsub_cols;
	
	for (i = 0; i < N_requested_rows; i++) {
		count_row = 0;
		subrow  = requested_rows[i];  /* row number in the local (sub) matrix */
		if (subrow < Nsub_rows)
			fullrow = sub_rows[subrow]; /* row number in the full matrix (of which
																		 this matrix is a part) */
		else {
			printf("getrow requested row %d of a submatrix with only %d rows\n",
						 subrow, Nsub_rows);
			exit(-1);
		}

		/* don't know the actual length of the returned row yet, but the length
			 of the row in the full matrix is a decent upper bound */
		row_lengths[i] = bindx[fullrow+1] - bindx[fullrow] + 1;
		if (count+row_lengths[i] > allocated_space) return(0);
		
		/* put in diagonal element if it exists */
		if (AZ_find_index(fullrow, sub_cols, Nsub_cols) >= 0) {
			columns[count + count_row  ] = subrow;
			values[count  + count_row++]  = val[fullrow];
		}


		/* off-diagonals */		
		for (j = bindx[fullrow] ; j < bindx[fullrow+1] ; j++) {
			fullcol = bindx[j];
			subcol = AZ_find_index(fullcol, sub_cols, Nsub_cols);
			if (subcol >= 0) {
				columns[count + count_row] = subcol;
				values[count  + count_row++] = val[j];
			}
		}
		/* now we know the real row length, so set it */
		row_lengths[i] = count_row;
		count += count_row;
	}
	return(1);
}

void AZ_subMSR_matvec_mult (double *b, double *c, struct AZ_MATRIX_STRUCT *Amat, 
														int proc_config[])
{
  double *val;
  int *data_org, *bindx, subrow, Nrows, Ncols, *rows, *cols;
  register int j, k, bindx_row;
  int          N, nzeros, fullrow, subcol;
	struct submat_struct *dataptr;


	dataptr = (struct submat_struct *) AZ_get_matvec_data(Amat);

	Nrows=dataptr->Nsub_rows;
	Ncols=dataptr->Nsub_cols;
	rows=dataptr->sub_rows;
	cols=dataptr->sub_cols;
  val = Amat->val;
  bindx = Amat->bindx;