1	/*//////////////////////////////////////////////////////////////////////////////
     2	// <file
     3	//   revision         = "$Revision: 1.1 $"
     4	//   date             = "$Date: 2009/02/17 16:45:17 $"
     5	//   tag              = "$Name: wn_sparse_solve-0-5 $"
     6	//   template_version = "cp_template.0.12"
     7	//   type             = "public"
     8	// >
     9	//   <license>
    10	//     See the src/README.txt file in this distribution for copyright and
    11	//     license information.
    12	//   </license>
    13	//   <description>
    14	//     <abstract>
    15	//       Example to demonstrate how to use wn_sparse_solve routines to 
    16	//       solve the linear matrix equation dY/dt = AY + P, where P is
    17	//       a constant vector, with Sparskit.
    18	//     </abstract>
    19	//   </description>
    20	// </file>
    21	//
    22	//////////////////////////////////////////////////////////////////////////////*/
    23	
    24	#include <WnSparseSolve.h>
    25	#include "solution_check.c"
    26	
    27	int
    28	main( int argc, char *argv[] )
    29	{
    30	
    31	  WnSparseSolve__Phi *p_my_phi;
    32	  WnMatrix *p_matrix; 
    33	  size_t i;
    34	  double d_t;
    35	  gsl_vector *p_initial_vector, *p_constant_vector, *p_my_solution_vector;
    36	
    37	  /*============================================================================
    38	  // Check input.
    39	  //==========================================================================*/
    40	
    41	  if( argc != 8 && argc != 9 ) {
    42	    fprintf(
    43	      stderr,
    44	      "\nUsage: %s matrix_file initial_file constant_vector_file time itmax workspace tol\n\n",
    45	      argv[0]
    46	    );
    47	    fprintf( stderr, "  matrix_file = name of input matrix file\n\n" );
    48	    fprintf(
    49	      stderr,
    50	      "  initial_file = name of input initial vector file\n\n"
    51	    );
    52	    fprintf( stderr,
    53	      "  constant_vector_file = name of input constant vector file\n\n"
    54	    );
    55	    fprintf( stderr, "  time = time over which to evolve\n\n" );
    56	    fprintf( stderr, "  itmax = maximum number of iterations\n\n" );
    57	    fprintf(
    58	      stderr,
    59	      "  workspace = integer size of workspace (0 < workspace < 50)\n\n"
    60	    );
    61	    fprintf( stderr, "  tol = tolerance (> 0)\n\n" );
    62	    fprintf( stderr, "  debug = debug (optional)\n\n" );
    63	    return EXIT_FAILURE;
    64	  }
    65	
    66	  /*============================================================================
    67	  // Assign time.
    68	  //==========================================================================*/
    69	
    70	  d_t = atof( argv[4] );
    71	
    72	  if( d_t < 0 ) {
    73	    printf( "Time should be > 0.\n" );
    74	    return EXIT_FAILURE;
    75	  }
    76	
    77	  /*============================================================================
    78	  // Get matrix.
    79	  //==========================================================================*/
    80	
    81	  p_matrix = WnMatrix__new_from_xml( argv[1], NULL );
    82	
    83	  /*============================================================================
    84	  // Get initial vector.
    85	  //==========================================================================*/
    86	
    87	  p_initial_vector =
    88	    WnMatrix__new_gsl_vector_from_xml( argv[2], NULL );
    89	
    90	  /*============================================================================
    91	  // Check that initial vector has the appropriate length.
    92	  //==========================================================================*/
    93	
    94	  if(
    95	    WnMatrix__get_gsl_vector_size( p_initial_vector )
    96	    != WnMatrix__getNumberOfColumns( p_matrix )
    97	  ) {
    98	    printf(
    99	      "\nInput initial vector does not have the right number of entries!\n"
   100	    );
   101	    return EXIT_FAILURE;
   102	  }
   103	
   104	  /*============================================================================
   105	  // Get the constant vector file.
   106	  //==========================================================================*/
   107	
   108	  p_constant_vector =
   109	    WnMatrix__new_gsl_vector_from_xml( argv[3], NULL );
   110	
   111	  /*============================================================================
   112	  // Check that constant vector has the appropriate length.
   113	  //==========================================================================*/
   114	
   115	  if(
   116	    WnMatrix__get_gsl_vector_size( p_constant_vector )
   117	    != WnMatrix__getNumberOfColumns( p_matrix )
   118	  ) {
   119	    printf(
   120	      "\nInput initial vector does not have the right number of entries!\n"
   121	    );
   122	    return EXIT_FAILURE;
   123	  }
   124	
   125	  /*============================================================================
   126	  // Get solver.
   127	  //==========================================================================*/
   128	
   129	  p_my_phi = WnSparseSolve__Phi__new( );
   130	
   131	  /*============================================================================
   132	  // Change settings.
   133	  //==========================================================================*/
   134	
   135	  WnSparseSolve__Phi__updateMaximumIterations( p_my_phi, atoi( argv[5] ) );
   136	
   137	  WnSparseSolve__Phi__updateWorkSpace( p_my_phi, atoi( argv[6] ) );
   138	
   139	  WnSparseSolve__Phi__updateTolerance( p_my_phi, atof( argv[7] ) );
   140	
   141	  /*============================================================================
   142	  // Set debug.
   143	  //==========================================================================*/
   144	
   145	  if( argv[8] ) {
   146	    if( strcmp( argv[8], "debug" ) == 0 )
   147	       WnSparseSolve__Phi__setDebug( p_my_phi );
   148	  }
   149	
   150	  /*============================================================================
   151	  // Get solution.
   152	  //==========================================================================*/
   153	
   154	  p_my_solution_vector =
   155	    WnSparseSolve__Phi__solve(
   156	      p_my_phi, p_matrix, p_initial_vector, p_constant_vector, d_t
   157	    );
   158	
   159	  /*============================================================================
   160	  // Check that the solution succeeded.
   161	  //==========================================================================*/
   162	
   163	  if( !is_valid_solution( p_my_solution_vector, argv[8] ) ) {
   164	    gsl_vector_free( p_initial_vector );
   165	    gsl_vector_free( p_constant_vector );
   166	    WnMatrix__free( p_matrix );
   167	    WnSparseSolve__Phi__free( p_my_phi );
   168	    return EXIT_FAILURE;
   169	  }
   170	
   171	  /*============================================================================
   172	  // Print out solution vector.
   173	  //==========================================================================*/
   174	
   175	  printf( "\nSolution vector:\n\n" );
   176	  for( i = 0; i < WnMatrix__getNumberOfRows( p_matrix ); i++ ) {
   177	    fprintf(
   178	      stdout,
   179	      "i = %lu  a_w[i] = %f\n",
   180	      (unsigned long) i,
   181	      gsl_vector_get( p_my_solution_vector, i )
   182	    );
   183	  }
   184	
   185	  /*============================================================================
   186	  // Clean up.
   187	  //==========================================================================*/
   188	 
   189	  gsl_vector_free( p_initial_vector );
   190	  gsl_vector_free( p_constant_vector );
   191	  gsl_vector_free( p_my_solution_vector );
   192	
   193	  WnMatrix__free( p_matrix );
   194	
   195	  WnSparseSolve__Phi__free( p_my_phi );
   196	
   197	  /*============================================================================
   198	  // Done.
   199	  //==========================================================================*/
   200	
   201	  return EXIT_SUCCESS;
   202	
   203	}