Ifpack Package Browser (Single Doxygen Collection) Development
distrib_msr_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_msr_matrix(int *proc_config, int *N_global,
00036         int *n_nonzeros, int *N_update, int **update,
00037         double **val, int **bindx,
00038         double **x, double **b, double **bt, double **xexact)
00039 #undef DEBUG 
00040 
00041 {
00042   int i, n_entries, N_columns, n_global_nonzeros;
00043   int ii, j, row, have_xexact = 0 ;
00044   int kk = 0;
00045   int max_ii = 0, max_jj = 0;
00046   int ione = 1;
00047   double value;
00048   double *cnt;
00049   int *pntr, *bindx1, *pntr1;
00050   double *val1, *b1, *bt1, *x1, *xexact1;
00051 
00052   printf("Processor %d of %d entering distrib_matrix.\n",
00053    proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00054 
00055   /*************** Distribute global matrix to all processors ************/
00056 
00057   if(proc_config[PAZ_node] == 0)
00058     {
00059       if ((*xexact) != NULL) have_xexact = 1;
00060       printf("Broadcasting exact solution\n");
00061     }
00062 
00063   if(proc_config[PAZ_N_procs]  > 1) { 
00064 
00065       PAZ_broadcast((char *) N_global,        sizeof(int), proc_config, PAZ_PACK);
00066       PAZ_broadcast((char *) n_nonzeros, sizeof(int), proc_config, PAZ_PACK);
00067       PAZ_broadcast((char *) &have_xexact, sizeof(int), proc_config, PAZ_PACK);
00068       PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
00069 
00070       if(proc_config[PAZ_node] != 0)
00071   {
00072     (*bindx) = (int   *) calloc(*n_nonzeros+1,sizeof(int)) ;
00073     (*val) = (double *) calloc(*n_nonzeros+1,sizeof(double)) ;
00074   }
00075 
00076       PAZ_broadcast((char *) (*bindx), sizeof(int)   *(*n_nonzeros+1), 
00077        proc_config, PAZ_PACK);
00078       PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
00079       PAZ_broadcast((char *) (*val),  sizeof(double)*(*n_nonzeros+1), 
00080        proc_config, PAZ_PACK);
00081       PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
00082 
00083       printf("Processor %d of %d done with matrix broadcast.\n",
00084        proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00085  
00086       /* Set rhs and initialize guess */
00087       if(proc_config[PAZ_node] != 0)
00088   {
00089     (*b) = (double *) calloc(*N_global,sizeof(double)) ;
00090     (*bt) = (double *) calloc(*N_global,sizeof(double)) ;
00091     (*x) = (double *) calloc(*N_global,sizeof(double)) ;
00092     if (have_xexact)
00093     (*xexact) =   (double *) calloc(*N_global,sizeof(double)) ;
00094   }
00095 
00096       PAZ_broadcast((char *) (*x), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
00097       PAZ_broadcast((char *) (*b), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
00098       PAZ_broadcast((char *) (*bt), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
00099       if (have_xexact)
00100   PAZ_broadcast((char *) 
00101          (*xexact), sizeof(double)*(*N_global), proc_config, PAZ_PACK);
00102       PAZ_broadcast(NULL, 0, proc_config, PAZ_SEND);
00103       printf("Processor %d of %d done with rhs/guess broadcast.\n",
00104        proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00105 
00106     }
00107 
00108   /********************** Generate update map  *************************/
00109 
00110   PAZ_read_update(N_update, update, proc_config, (*N_global),
00111            1, PAZ_linear) ;
00112   
00113   printf("Processor %d of %d has %d rows of %d total rows.\n",
00114    proc_config[PAZ_node],proc_config[PAZ_N_procs],*N_update,(*N_global)) ;
00115 
00116   /*************** Construct local matrix from global matrix ************/
00117 
00118   /* The local matrix is a copy of the rows assigned to this processor.  
00119      It is stored in MSR format and still has global indices (PAZ_transform
00120      will complete conversion to local indices.
00121   */
00122 
00123   if(proc_config[PAZ_N_procs]  > 1) { 
00124       n_global_nonzeros = *n_nonzeros;
00125 
00126       *n_nonzeros = *N_update;
00127       
00128       for (i=0; i<*N_update; i++)
00129   *n_nonzeros += (*bindx)[(*update)[i]+1] - (*bindx)[(*update)[i]];
00130 
00131       printf("Processor %d of %d has %d nonzeros of %d total nonzeros.\n",
00132        proc_config[PAZ_node],proc_config[PAZ_N_procs],
00133        *n_nonzeros,n_global_nonzeros) ;
00134 
00135 #ifdef DEBUG
00136       { double sum1 = 0.0;
00137       for (i=0;i<(*N_global); i++) sum1 += (*b)[i];
00138 
00139       printf("Processor %d of %d has sum of b = %12.4g.\n",
00140        proc_config[PAZ_node],proc_config[PAZ_N_procs],sum1) ;
00141       }
00142 #endif /* DEBUG */
00143 
00144       /* Allocate memory for local matrix */
00145 
00146       bindx1 = (int   *) calloc(*n_nonzeros+1,sizeof(int)) ;
00147       val1 = (double *) calloc(*n_nonzeros+1,sizeof(double)) ;
00148       b1 =   (double *) calloc(*N_update,sizeof(double)) ;
00149       bt1 =   (double *) calloc(*N_update,sizeof(double)) ;
00150       x1 =   (double *) calloc(*N_update,sizeof(double)) ;
00151       if (have_xexact)
00152       xexact1 =   (double *) calloc(*N_update,sizeof(double)) ;
00153      
00154       bindx1[0] = *N_update+1;
00155       
00156       for (i=0; i<*N_update; i++)
00157   {
00158     row = (*update)[i];
00159     b1[i] = (*b)[row];
00160     bt1[i] = (*bt)[row];
00161     x1[i] = (*x)[row];
00162     if (have_xexact) xexact1[i] = (*xexact)[row];
00163     val1[i] = (*val)[row];
00164     bindx1[i+1] = bindx1[i];
00165 
00166 #ifdef DEBUG    
00167     printf("Proc %d of %d: Global row = %d: Local row = %d: 
00168                   b = %12.4g: x = %12.4g: bindx = %d: val = %12.4g \n",
00169      proc_config[PAZ_node],proc_config[PAZ_N_procs], 
00170      row, i, b1[i], x1[i], bindx1[i], val1[i]) ;
00171 #endif
00172 
00173     for (j = (*bindx)[row]; j < (*bindx)[row+1]; j++)
00174       {
00175         val1[  bindx1 [i+1] ] = (*val)[j];
00176         bindx1[bindx1 [i+1] ] = (*bindx)[j];
00177         bindx1[i+1] ++;
00178       }
00179   }
00180 
00181       printf("Processor %d of %d done with extracting local operators.\n",
00182        proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00183 
00184       if (have_xexact)
00185   {
00186     printf(
00187      "The residual using MSR format and exact solution on processor %d is %12.4g\n",
00188         proc_config[PAZ_node],
00189         smsrres (*N_update, (*N_global), val1, bindx1, xexact1, (*xexact), b1));
00190   }
00191   
00192       /* Release memory for global matrix, rhs and solution */
00193       
00194       free ((void *) (*val));
00195       free ((void *) (*bindx));
00196       free ((void *) (*b));
00197       free ((void *) (*bt));
00198       free ((void *) (*x));
00199       if (have_xexact) free((void *) *xexact);
00200 
00201       /* Return local matrix through same pointers. */
00202       
00203       *val = val1;
00204       *bindx = bindx1;
00205       *b = b1;
00206       *bt = bt1;
00207       *x = x1;
00208       if (have_xexact) *xexact = xexact1;
00209 
00210   }
00211   if (have_xexact && proc_config[PAZ_N_procs]  == 1)
00212     {
00213       printf(
00214        "The residual using MSR format and exact solution on processor %d is %12.4g\n",
00215        proc_config[PAZ_node],
00216        smsrres (*N_update, (*N_global), (*val), (*bindx), 
00217           (*xexact), (*xexact), (*b)));
00218     }
00219   
00220   
00221   printf("Processor %d of %d leaving distrib_matrix.\n",
00222    proc_config[PAZ_node],proc_config[PAZ_N_procs]) ;
00223   
00224   /* end distrib_matrix */
00225 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines