Ifpack Package Browser (Single Doxygen Collection) Development
distrib_vbr_matrix.c
Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include <stdlib.h>
00031 #include <stdio.h>
00032 #include "paz_aztec.h"
00033 #include "prototypes.h"
00034 
00035 void distrib_vbr_matrix(int *proc_config,
00036         int *N_global, int *N_blk_global,
00037         int *n_nonzeros, int *n_blk_nonzeros, 
00038         int *N_update, int **update,
00039         double **val, int **indx, int **rpntr, int **cpntr,
00040         int **bpntr, int **bindx,
00041         double **x, double **b, double **bt, double **xexact)
00042 #undef DEBUG 
00043 
00044 {
00045   int i, n_entries, N_columns, n_global_nonzeros, n_global_blk_nonzeros;
00046   int N_local;
00047   int ii, j, row, have_xexact = 0 ;
00048   int kk = 0;
00049   int max_ii = 0, max_jj = 0;
00050   int ione = 1;
00051   double value;
00052   double *cnt;
00053   int *rpntr1, *bindx1, *bpntr1, *indx1;
00054   double *val1, *b1, *bt1, *x1, *xexact1;
00055 
00056   printf("Processor %d of %d entering distrib_matrix.\n",
00057    proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00058 
00059   /*************** Distribute global matrix to all processors ************/
00060 
00061   if(proc_config[PAZ_node] == 0)
00062     {
00063       if ((*xexact) != NULL) have_xexact = 1;
00064       printf("Broadcasting exact solution\n");
00065     }
00066 
00067   if(proc_config[PAZ_N_procs]  > 1)
00068     { 
00069 
00070       PAZ_broadcast((char *) N_global,      sizeof(int), proc_config, PAZ_PACK);
00071       PAZ_broadcast((char *) N_blk_global,  sizeof(int), proc_config, PAZ_PACK);
00072       PAZ_broadcast((char *) n_nonzeros,     sizeof(int), proc_config, PAZ_PACK);
00073       PAZ_broadcast((char *) n_blk_nonzeros, sizeof(int), proc_config, PAZ_PACK);
00074       PAZ_broadcast((char *) &have_xexact,   sizeof(int), proc_config, PAZ_PACK);
00075       PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
00076 
00077       printf("Processor %d of %d done with global parameter  broadcast.\n",
00078        proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00079 
00080       if(proc_config[PAZ_node] != 0)
00081   {
00082       *bpntr = (int   *) calloc(*N_blk_global+1,sizeof(int)) ;
00083       *rpntr = (int   *) calloc(*N_blk_global+1,sizeof(int)) ;
00084       *bindx = (int   *) calloc(*n_blk_nonzeros+1,sizeof(int)) ;
00085       *indx  = (int   *) calloc(*n_blk_nonzeros+1,sizeof(int)) ;
00086       *val = (double *) calloc(*n_nonzeros+1,sizeof(double)) ;
00087       printf("Processor %d of %d done with global calloc.\n",
00088        proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00089 }
00090 
00091       PAZ_broadcast((char *) (*bpntr), sizeof(int)   *(*N_blk_global+1), 
00092        proc_config, PAZ_PACK);
00093       PAZ_broadcast((char *) (*rpntr), sizeof(int)   *(*N_blk_global+1), 
00094        proc_config, PAZ_PACK);
00095       PAZ_broadcast((char *) (*bindx), sizeof(int)   *(*n_blk_nonzeros+1), 
00096        proc_config, PAZ_PACK);
00097       PAZ_broadcast((char *) (*indx),  sizeof(int)   *(*n_blk_nonzeros+1), 
00098        proc_config, PAZ_PACK);
00099       PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
00100       PAZ_broadcast((char *) (*val),  sizeof(double)*(*n_nonzeros+1), 
00101        proc_config, PAZ_PACK);
00102       PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
00103 
00104       printf("Processor %d of %d done with matrix broadcast.\n",
00105        proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00106  
00107       /* Set rhs and initialize guess */
00108       if(proc_config[PAZ_node] != 0)
00109   {
00110     (*b) = (double *) calloc(*N_global,sizeof(double)) ;
00111     (*bt) = (double *) calloc(*N_global,sizeof(double)) ;
00112     (*x) = (double *) calloc(*N_global,sizeof(double)) ;
00113     if (have_xexact)
00114     (*xexact) =   (double *) calloc(*N_global,sizeof(double)) ;
00115   }
00116 
00117       PAZ_broadcast((char *) (*x), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
00118       PAZ_broadcast((char *) (*b), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
00119       PAZ_broadcast((char *) (*bt), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
00120       if (have_xexact)
00121   PAZ_broadcast((char *) 
00122          (*xexact), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
00123       PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
00124       printf("Processor %d of %d done with rhs/guess broadcast.\n",
00125        proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00126 
00127     }
00128 
00129   /********************** Generate update map  *************************/
00130 
00131   PAZ_read_update(N_update, update, proc_config, *N_blk_global,
00132            1, PAZ_linear) ;
00133 
00134   printf("Processor %d of %d has %d rows of %d total block rows.\n",
00135    proc_config[PAZ_node],proc_config[PAZ_N_procs],*N_update,*N_blk_global) ;
00136 
00137   /*************** Construct local matrix from global matrix ************/
00138 
00139   /* The local matrix is a copy of the rows assigned to this processor.  
00140      It is stored in MSR format and still has global indices (PAZ_transform
00141      will complete conversion to local indices.
00142   */
00143 
00144   if(proc_config[PAZ_N_procs]  > 1)
00145     { 
00146       n_global_nonzeros = *n_nonzeros;
00147       n_global_blk_nonzeros = *n_blk_nonzeros;
00148 
00149       *n_nonzeros = 0;
00150       *n_blk_nonzeros = 0;
00151       N_local = 0;
00152       
00153       for (i=0; i<*N_update; i++)
00154   {
00155     row = (*update)[i];
00156     *n_nonzeros     += (*indx)[(*bpntr)[row+1]] - (*indx)[(*bpntr)[row]];
00157     *n_blk_nonzeros += (*bpntr)[row+1] - (*bpntr)[row];
00158     N_local         += (*rpntr)[row+1] - (*rpntr)[row];
00159     
00160   }
00161 
00162       printf("Processor %d of %d has %d nonzeros of %d total nonzeros.\n",
00163        proc_config[PAZ_node],proc_config[PAZ_N_procs],
00164        *n_nonzeros,n_global_nonzeros) ;
00165 
00166    printf("Processor %d of %d has %d block nonzeros of %d total block nonzeros.\n",
00167        proc_config[PAZ_node],proc_config[PAZ_N_procs],
00168        *n_blk_nonzeros,n_global_blk_nonzeros) ;
00169 
00170    printf("Processor %d of %d has %d equations of %d total equations.\n",
00171        proc_config[PAZ_node],proc_config[PAZ_N_procs],
00172        N_local,*N_global) ;
00173 
00174 #ifdef DEBUG
00175       { double sum1 = 0.0;
00176       for (i=0;i<*N_global; i++) sum1 += (*b)[i];
00177 
00178       printf("Processor %d of %d has sum of b = %12.4g.\n",
00179        proc_config[PAZ_node],proc_config[PAZ_N_procs],sum1) ;
00180       }
00181 #endif /* DEBUG */
00182 
00183       /* Allocate memory for local matrix */
00184 
00185       bpntr1 = (int   *) calloc(*N_update+1,sizeof(int)) ;
00186       rpntr1 = (int   *) calloc(*N_update+1,sizeof(int)) ;
00187       bindx1 = (int   *) calloc(*n_blk_nonzeros+1,sizeof(int)) ;
00188       indx1  = (int   *) calloc(*n_blk_nonzeros+1,sizeof(int)) ;
00189       val1 = (double *) calloc(*n_nonzeros+1,sizeof(double)) ;
00190       b1 =   (double *) calloc(N_local,sizeof(double)) ;
00191       bt1 =   (double *) calloc(N_local,sizeof(double)) ;
00192       x1 =   (double *) calloc(N_local,sizeof(double)) ;
00193       if (have_xexact)
00194       xexact1 =   (double *) calloc(N_local,sizeof(double)) ;
00195 
00196       {     
00197   int cur_blk_size, indx_offset, len_val, row_offset, row_offset1;
00198   double *val_ptr, *val1_ptr;
00199 
00200   bpntr1[0] = 0;
00201   indx1[0] = 0;
00202   rpntr1[0] = 0;
00203   for (i=0; i<*N_update; i++)
00204     {
00205       row = (*update)[i];
00206       cur_blk_size = (*rpntr)[row+1] - (*rpntr)[row];
00207       rpntr1[i+1] = rpntr1[i] + cur_blk_size;
00208       row_offset = (*rpntr)[row];
00209       row_offset1 = rpntr1[i];
00210       for (j = 0; j<cur_blk_size; j++)
00211         {
00212     b1[row_offset1+j] = (*b)[row_offset+j];
00213     x1[row_offset1+j] = (*x)[row_offset+j];
00214     if (have_xexact) xexact1[row_offset1+j] = (*xexact)[row_offset+j];
00215         }
00216       bpntr1[i+1] = bpntr1[i];
00217       
00218 #ifdef DEBUG    
00219       printf("Proc %d of %d: Global row = %d: Local row = %d: 
00220                     b = %12.4g: x = %12.4g: bindx = %d: val = %12.4g \n",
00221         proc_config[PAZ_node],proc_config[PAZ_N_procs], 
00222         row, i, b1[i], x1[i], bindx1[i], val1[i]) ;
00223 #endif
00224       indx_offset = (*indx)[(*bpntr)[row]] - indx1[bpntr1[i]];
00225       for (j = (*bpntr)[row]; j < (*bpntr)[row+1]; j++)
00226         {
00227     indx1[bpntr1 [i+1] + 1] = (*indx)[j+1] - indx_offset;
00228     bindx1[bpntr1 [i+1] ] = (*bindx)[j];
00229     bpntr1[i+1] ++;
00230         }
00231       len_val = indx1[bpntr1[i+1]] - indx1[bpntr1[i]];
00232       val_ptr = (*val)+(*indx)[(*bpntr)[row]];
00233       val1_ptr = val1+indx1[bpntr1[i]];
00234       for (j = 0; j<len_val; j++)
00235         { 
00236     *val1_ptr = *val_ptr;
00237     val_ptr++; val1_ptr++;
00238         }
00239     }
00240       }
00241       printf("Processor %d of %d done with extracting local operators.\n",
00242        proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00243 
00244       if (have_xexact)
00245   {
00246     printf(
00247      "The residual using VBR format and exact solution on processor %d is %12.4g\n",
00248         proc_config[PAZ_node],
00249         svbrres (N_local, *N_global, *N_update, val1, indx1, bindx1, 
00250            rpntr1, (*rpntr), bpntr1, bpntr1+1,
00251            (*xexact), b1));
00252   }
00253   
00254       /* Release memory for global matrix, rhs and solution */
00255       
00256       free ((void *) (*val));
00257       free ((void *) (*indx));
00258       free ((void *) (*bindx));
00259       free ((void *) (*bpntr));
00260       free ((void *) (*rpntr));
00261       free ((void *) (*b));
00262       free ((void *) (*bt));
00263       free ((void *) (*x));
00264       if (have_xexact) free((void *) *xexact);
00265 
00266       /* Return local matrix through same pointers. */
00267       
00268       *val = val1;
00269       *indx = indx1;
00270       *bindx = bindx1;
00271       *bpntr = bpntr1;
00272       *rpntr = rpntr1;
00273       *b = b1;
00274       *bt = bt1;
00275       *x = x1;
00276       if (have_xexact) *xexact = xexact1;
00277 
00278     }
00279       if (have_xexact && proc_config[PAZ_N_procs]  == 1)
00280   {
00281     printf(
00282      "The residual using VBR format and exact solution on processor %d is %12.4g\n",
00283         proc_config[PAZ_node],
00284         svbrres (*N_global, *N_global, *N_update, (*val), (*indx), (*bindx), 
00285            (*rpntr), (*rpntr), (*bpntr), (*bpntr)+1,
00286            (*xexact), (*b)));
00287   }
00288 
00289   
00290   printf("Processor %d of %d leaving distrib_matrix.\n",
00291    proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00292   
00293   /* end distrib_matrix */
00294 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines