az_matrix_util.c

并行解法器,功能强大

  data_org = Amat->data_org;


  N = data_org[AZ_N_internal] + data_org[AZ_N_border];

  /* exchange boundary info */

  AZ_exchange_bdry(b, data_org, proc_config);

  for (subrow = 0; subrow < Nrows; subrow++) {
		fullrow = rows[subrow];  /* this row in the full matrix */

    /* compute diagonal contribution if there is one*/
		if (AZ_find_index(fullrow, cols, Ncols) >= 0)
			*c = val[fullrow] * b[subrow];
		else
			*c = 0.0;

    /* nonzero off diagonal contribution */

    bindx_row = bindx[fullrow];
    nzeros    = bindx[fullrow+1] - bindx_row;

		/* loop through the nonzeros in the full matrix and check to see if they're
			 in the submatrix.  If they are, add the appropriate product to the appropriate
			 entry */
    for (j = 0; j < nzeros; j++) {
      k   = bindx_row + j;
			subcol = AZ_find_index(bindx[k], cols, Ncols);
			if (subcol >= 0)				
				*c += val[k] * b[subcol];
    }
    c++;
  }
 
} /* AZ_subMSR_matvec_mult */



struct blockmat_struct {
	int Nblock_rows, Nblock_cols;
	int *Nsub_rows, **sub_rows, *Nsub_cols, **sub_cols, Nsub_mats;
	AZ_MATRIX **submat_list;
	int **submat_locs;
	int tot_Nrows;
};


AZ_MATRIX *AZ_blockmatrix_create(AZ_MATRIX **submat_list, int Nsub_mats, int **submat_locs, 
				int Nblock_rows, int Nblock_cols, int Nsub_rows[], int **sub_rows, int Nsub_cols[], 
				int **sub_cols, int proc_config[])
		 /* submat_locs is a Nsub_mats x 2 array of integers - the first column gives the block
				row of each submatrix and the second column gives the block column of each submatrix */
{
	int i, j, tot_Nrows;
	
	AZ_MATRIX *block_matrix;
	struct blockmat_struct *new_dataptr;

	
	tot_Nrows=0;
	for (i = 0; i < Nblock_rows; i++)
		tot_Nrows += Nsub_rows[i];
	
	block_matrix = AZ_matrix_create(tot_Nrows);
	
	new_dataptr = (struct blockmat_struct *)malloc(sizeof(struct blockmat_struct));
	
	new_dataptr->tot_Nrows   = tot_Nrows;
	new_dataptr->Nblock_rows = Nblock_rows;
	new_dataptr->Nblock_cols = Nblock_cols;
	new_dataptr->Nsub_mats   = Nsub_mats;

	/* allocate space for all the things we need to copy over */
	new_dataptr->submat_list = (AZ_MATRIX **) malloc(Nsub_mats * sizeof (AZ_MATRIX *));
	new_dataptr->submat_locs = (int **)  malloc(Nsub_mats   * sizeof(int *));
	new_dataptr->Nsub_rows   = (int *)  malloc(Nblock_rows * sizeof(int));
	new_dataptr->Nsub_cols   = (int *)  malloc(Nblock_cols * sizeof(int));
	new_dataptr->sub_rows    = (int **) malloc(Nblock_rows * sizeof(int *));
	new_dataptr->sub_cols    = (int **) malloc(Nblock_cols * sizeof(int *));
	
	if (new_dataptr->sub_cols == NULL) {
		printf("memory allocation error\n");
		exit(-1);
	}
	
	for (i = 0; i < Nsub_mats; i++) {
		new_dataptr->submat_list[i]=submat_list[i];
		new_dataptr->submat_locs[i]=(int *) malloc(2*sizeof(int));
		if (new_dataptr->submat_locs[i] == NULL) {
			printf("memory allocation error\n");
			exit(-1);
		}
		new_dataptr->submat_locs[i][0]=submat_locs[i][0];
		new_dataptr->submat_locs[i][1]=submat_locs[i][1];
	}

	for (i = 0; i < Nblock_rows; i++) {
		new_dataptr->Nsub_rows[i] = Nsub_rows[i];
		new_dataptr->sub_rows[i] = (int *) malloc(Nsub_rows[i] * sizeof(int));
		if (new_dataptr->sub_rows[i] == NULL) {
			printf("memory allocation error\n");
			exit(-1);
		}
		
		for (j = 0; j < Nsub_rows[i]; j++)
			new_dataptr->sub_rows[i][j]=sub_rows[i][j];
	}

	for (i = 0; i < Nblock_cols; i++) {
		new_dataptr->Nsub_cols[i] = Nsub_cols[i];
		new_dataptr->sub_cols[i] = (int *) malloc(Nsub_cols[i] * sizeof(int));
		if (new_dataptr->sub_cols[i] == NULL) {
			printf("memory allocation error\n");
			exit(-1);
		}
		
		for (j = 0; j < Nsub_cols[i]; j++)
			new_dataptr->sub_cols[i][j]=sub_cols[i][j];
	}
			
	
	AZ_set_MATFREE(block_matrix, new_dataptr, AZ_blockMSR_matvec_mult);
	AZ_set_MATFREE_getrow(block_matrix, new_dataptr, AZ_blockMSR_getrow, NULL, 0, proc_config);


	return(block_matrix);
}

void AZ_blockmatrix_destroy(AZ_MATRIX **blockmat)
{
	int **sub_rows, **sub_cols, *Nsub_rows, *Nsub_cols, **submat_locs;
	int i, Nblock_rows, Nblock_cols, Nsub_mats;
	struct blockmat_struct *data;
	AZ_MATRIX **submat_list;

	data = (struct blockmat_struct *) AZ_get_matvec_data(*blockmat);

	sub_rows = data->sub_rows;
	sub_cols = data->sub_cols;
	Nsub_rows = data->Nsub_rows;
	Nsub_cols = data->Nsub_cols;
	Nblock_rows = data->Nblock_rows;
	Nblock_cols = data->Nblock_cols;
	Nsub_mats = data->Nsub_mats;
	submat_locs = data->submat_locs;
	submat_list = data->submat_list;

	for (i=0; i < Nblock_rows; i++)
		free(sub_rows[i]);
	for (i=0; i < Nblock_cols; i++)
		free(sub_cols[i]);
	for (i=0; i < Nsub_mats; i++)
		free(submat_locs[i]);

	free(Nsub_rows);
	free(Nsub_cols);
	free(submat_list);
	free(submat_locs);
	free(sub_rows);
	free(sub_cols);
	free(data);
	AZ_matrix_destroy(blockmat);
}


int  AZ_blockMSR_getrow(int cols[], double vals[], int row_lengths[],
											AZ_MATRIX *Amat, int N_requested_rows,
											int requested_rows[], int allocated_space)
{ 
	int    i, j, k, ctr, count = 0, block_row, block_col, count_row;
	int     *Nsub_rows, *Nsub_cols, **sub_rows, **sub_cols, full_req_row;
	int    local_req_row, tmp_row_len, *tmpcols, Nsub_mats, **submat_locs;
	int    tmp_alloc_space, max_nnz_per_row=500, tmp;
	double *tmpvals;
	struct blockmat_struct *dataptr;
	struct AZ_MATRIX_STRUCT *submat;

	dataptr = (struct blockmat_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;
	Nsub_mats=dataptr->Nsub_mats;
	submat_locs=dataptr->submat_locs;

	tmpcols = (int *) malloc(max_nnz_per_row*sizeof(int));
	tmpvals = (double *) malloc(max_nnz_per_row*sizeof(double));
	if (tmpvals==NULL) {
		printf("memory allocation error\n");
		exit(-1);
	}
	tmp_alloc_space = max_nnz_per_row;


	for (i = 0; i < N_requested_rows; i++) {
		count_row = 0;
		full_req_row  = requested_rows[i];
		if (full_req_row > dataptr->tot_Nrows) {
			printf("Error: requested row %d of a matrix with %d rows\n", 
						 full_req_row+1, dataptr->tot_Nrows);
			exit(-1);
		}
			
		ctr=0;
		/* figure out which submatrix this row is in and which row it is 
			 within that submatrix */
		local_req_row=AZ_find_index(full_req_row, sub_rows[0],Nsub_rows[0]);
		while(local_req_row < 0) {
			ctr ++;
			local_req_row = AZ_find_index(full_req_row, sub_rows[ctr],
																		Nsub_rows[ctr]);			
		}

		block_row=ctr;  /* the block row where the requesed row is located */

		/* now loop through the submatrices and see which ones are in this block row: */
		for (j = 0; j < Nsub_mats; j++)
			if (submat_locs[j][0]==block_row) { /* this matrix contains part of the requested row */
				submat = dataptr->submat_list[j];
				/* figure out which block column this submatrix is in and get the relavent row */
				block_col=submat_locs[j][1];
				tmp = submat->getrow(tmpcols, tmpvals, &tmp_row_len, submat, 1, 
														 &local_req_row, tmp_alloc_space);
				/* in case we didn't allocate enough space for this row, keep checking the return
					 value of getrow until we get a positive result */
				while (tmp == 0) {
					free(tmpcols);
					free(tmpvals);
					max_nnz_per_row=max_nnz_per_row*2+1;
					tmp_alloc_space=max_nnz_per_row;
					tmpcols=(int *)malloc(tmp_alloc_space*sizeof(int));
					tmpvals=(double *)malloc(tmp_alloc_space*sizeof(double));
					tmp = submat->getrow(tmpcols, tmpvals, &tmp_row_len, submat, 1, 
															 &local_req_row, tmp_alloc_space);
				}

				for (k = 0; k < tmp_row_len; k++) {
					/* if we still have space in the output vector, go through and insert each
						 value from the submatrix row into the output vector */
					if (allocated_space > count+count_row) {
						cols[count+count_row]   = sub_cols[block_col][tmpcols[k]];
						vals[count+count_row++] = tmpvals[k];
					}
					else {
						free(tmpcols);
						free(tmpvals);
						return(0);
					}
				}
			}
		count += count_row;
		row_lengths[i]=count_row;
	}  

	free(tmpcols);
	free(tmpvals);
	return(1);
}

void AZ_blockMSR_matvec_mult (double *b, double *c, struct AZ_MATRIX_STRUCT *Amat, 
														int proc_config[])
{
	double *tmpb, *tmpc;
  int *data_org, Nrows;
  register int j;
  int          i, N;
	int *submat_loc, block_row, block_col, Nsub_mats, Nsub_rows, Nsub_cols;
	struct blockmat_struct *dataptr;
	AZ_MATRIX *submat;

  data_org = Amat->data_org;

  N = data_org[AZ_N_internal] + data_org[AZ_N_border];

  /* exchange boundary info */

  AZ_exchange_bdry(b, data_org, proc_config);


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

	Nrows=dataptr->tot_Nrows;
	tmpb= (double *) malloc(Nrows * sizeof(double));
	tmpc= (double *) malloc(Nrows * sizeof(double));
	if (tmpc == NULL) {
		printf("memory allocation error\n");
		exit(-1);
	}

	for (i = 0; i < Nrows; i++)
		c[i] = 0.0;
	
	Nsub_mats= dataptr->Nsub_mats;

	/* loop through all the submatrices */
	for (i = 0; i < Nsub_mats; i++) {
		submat = dataptr->submat_list[i];
		submat_loc=dataptr->submat_locs[i];
		block_row = submat_loc[0];
		block_col = submat_loc[1];
		Nsub_rows = dataptr->Nsub_rows[block_row];
		Nsub_cols = dataptr->Nsub_cols[block_col];
		/* copy the rows of b that correspond to the columns of this
			 submatrix into a temp vector so we can run matvec on it */
		for (j = 0; j < Nsub_cols; j++)
			tmpb[j]=b[dataptr->sub_cols[block_col][j]];

		submat->matvec(tmpb, tmpc, submat, proc_config);

		/* now add the result to the appropriate rows of the output
			 vector */
		for (j = 0; j < Nsub_rows; j++)
			c[dataptr->sub_rows[block_row][j]] += tmpc[j];
	}
 
} /* AZ_blockMSR_matvec_mult */


void AZ_abs_matvec_mult(double *b, double *c,AZ_MATRIX *Amat,int proc_config[])


/******************************************************************************
  c = |A||b|:
  Sparse (square) overlapped matrix-vector multiply, using the  MSR
  data structure .

  Author:          Lydie Prevost, SNL, 1422
  =======

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

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


  b:               Contains the vector b.

  c:               Contains the result vector c.



  options:         Determines specific solution method and other parameters.

  params:          Drop tolerance and convergence tolerance info.

  Amat:            Structure used to represent the matrix (see file az_aztec.h
                   and Aztec User's Guide).

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

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



{
  double *val, *x;
  int *data_org, *bindx;
  register int j, k, irow, bindx_row;
  int          N, nzeros, num_blks, bpoff, rpoff, iblk_row, ibpntr, 
               ibpntr_next = 0;
  int          m1, n1, irpntr, irpntr_next, ib1, ii, jj, iblk_size, jblk;
  double       *val_pntr, *c_pntr;
  int          *rpntr, *cpntr, *bpntr;

  val      = Amat->val; 
  bindx    = Amat->bindx;
  data_org = Amat->data_org;
  N = data_org[AZ_N_internal] + data_org[AZ_N_border];

  /* exchange boundary info */


  AZ_exchange_bdry(b, data_org, proc_config);

  if (data_org[AZ_matrix_type] == AZ_MSR_MATRIX) {
     for (irow = 0; irow < N; irow++) {
        *c = fabs(val[irow]) * fabs(b[irow]);
        bindx_row = bindx[irow];
        nzeros    = bindx[irow+1] - bindx_row;
        for (j = 0; j < nzeros; j++) {
           k   = bindx_row + j;
           *c += fabs(val[k]) * fabs(b[bindx[k]]);
        }
        c++;
     }
  }
  else if (data_org[AZ_matrix_type] == AZ_VBR_MATRIX) {
     cpntr = Amat->cpntr;
     rpntr = Amat->rpntr;
     bpntr = Amat->bpntr;
     num_blks = data_org[AZ_N_int_blk] + data_org[AZ_N_bord_blk];

     bpoff    = *bpntr;
     rpoff    = *rpntr;
     val_pntr = val;
     for (j = 0; j < rpntr[num_blks] - rpoff; c[j++] = 0.0);
     irpntr_next = *rpntr++;
     bpntr++;
     c          -= rpoff;

     for (iblk_row = 0; iblk_row < num_blks; iblk_row++) {
        irpntr      = irpntr_next;
        irpntr_next = *rpntr++;
        ibpntr      = ibpntr_next;
        ibpntr_next = *bpntr++ - bpoff;
        c_pntr      = c + irpntr;
        m1          = irpntr_next - irpntr;

        /* loop over each block in the current row block */

        for (j = ibpntr; j < ibpntr_next; j++) {
           jblk = *(bindx+j);
           ib1  = *(cpntr+jblk);
           n1   = cpntr[jblk + 1] - ib1;
           iblk_size = m1*n1;
           x    = b + ib1;
           for (ii = 0; ii < m1; ii++) {
              for (jj = 0; jj < n1; jj++) {
                 c_pntr[ii] += fabs(val_pntr[n1*jj + ii])*fabs(x[jj]);
              }
           }
           val_pntr += iblk_size;
       }
     }


  }
  else {
     printf("Error: AZ_expected_value convergence options can only be done with MSR or VBR matrices\n");
     AZ_exit(1);
  }
 
} /* AZ_abs_matvec_mult */