/*************************************************************************** Parallel Linear Algebra Package Release R2.0.2 6 Feb 2000 Copyright (c) 1997,1998,1999,2000 Robert van de Geijn and The University of Texas at Austin. See the file README for details on the gnu license This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 1, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. Written under the direction of: Robert van de Geijn, Department of Computer Sciences, University of Texas at Austin 78712. rvdg@cs.utexas.edu Many people contributed to the Parallel Linear Algebra Package (PLAPACK) including: Gregory A. Baker, Center for Space Research, University of Texas at Austin 78722 baker@csr.utexas.edu Philip A. Alpatov, Department of Physics, University of Texas at Austin 78712. philip@physics.utexas.edu James Overfelt, Department of Mathematics, University of Texas at Austin 78712. overfelt@math.utexas.edu John Andrew Gunnels, Department of Computer Sciences, University of Texas at Austin 78712. gunnels@cs.utexas.edu Greg Morrow, Department of Physics, University of Texas at Austin 78712. morrow@physics.utexas.edu Wesley Reiley, Department of Computer Sciences University of Texas at Austin 78712. wesley@cs.utexas.edu Dr. Robert Van de Geijn, Department of Computer Sciences, University of Texas at Austin 78712. rvdg@cs.utexas.edu This release (Release 2.0.2) was written primarily by Robert van de Geijn ***********i****************************************************************/ #include "PLA.h" static PLA_Obj alpha_cpy = NULL, A_cpy = NULL, x_cpy = NULL, beta_cpy = NULL, y_cpy = NULL; static int old_size_malloced; int PLA_Gemv_enter( int transa, PLA_Obj alpha, PLA_Obj A, PLA_Obj x, PLA_Obj beta, PLA_Obj y ) { int value = PLA_SUCCESS, size, length_A, width_A, length_x, width_x, length_y, width_y, objtype, proj_onto; char routine_name[ 35 ] = "PLA_Gemv"; PLA_Routine_stack_push( routine_name ); PLA_Routine_stack_push( "PLA_Gemv_enter" ); old_size_malloced = PLA_Total_size_malloced( ); if ( PLA_CHECK_PARAMETERS ){ /* Check if transpose parameters are valid */ if ( !PLA_Valid_trans_parameter( transa ) ){ PLA_Warning( "Invalid parameter transa" ); value--; } /* Check if alpha is valid multiscalar of size 1x1 */ if ( alpha == NULL || !PLA_Valid_object( alpha ) ) { PLA_Warning( "Invalid object alpha" ); value--; } PLA_Obj_objtype( alpha, &objtype ); if ( objtype != PLA_MSCALAR ){ PLA_Warning( "Invalid objtype for alpha" ); value--; } PLA_Obj_global_length( alpha, &size ); if ( size != 1 ){ PLA_Warning( "Invalid global length for alpha" ); value--; } PLA_Obj_global_width( alpha, &size ); if ( size != 1 ){ PLA_Warning( "Invalid global width for alpha" ); value--; } /* Check if A is valid matrix */ if ( A == NULL || !PLA_Valid_object( A ) ) { PLA_Warning( "Invalid object A" ); value--; } PLA_Obj_objtype( A, &objtype ); if ( objtype != PLA_MATRIX ){ PLA_Warning( "Invalid objtype for A" ); value--; } PLA_Obj_global_length( A, &length_A ); PLA_Obj_global_width ( A, &width_A ); /* Check if x is valid vector */ if ( x == NULL || !PLA_Valid_object( x ) ) { PLA_Warning( "Invalid object x" ); value--; } PLA_Obj_objtype( x, &objtype ); if ( objtype != PLA_MATRIX && objtype != PLA_MVECTOR && objtype != PLA_PMVECTOR ){ PLA_Warning( "Invalid objtype for x" ); value--; } PLA_Obj_project_onto( x, &proj_onto ); if ( proj_onto == PLA_PROJ_ONTO_COL ){ PLA_Obj_global_length( x, &length_x ); PLA_Obj_global_width( x, &width_x ); } else { PLA_Obj_global_width( x, &length_x ); PLA_Obj_global_length( x, &width_x ); } if ( width_x != 1 ){ PLA_Warning( "x is not of width 1" ); value--; } /* Check if beta is valid multiscalar of size 1x1 */ if ( beta == NULL || !PLA_Valid_object( beta ) ) { PLA_Warning( "Invalid object beta" ); value--; } PLA_Obj_objtype( beta, &objtype ); if ( objtype != PLA_MSCALAR ){ PLA_Warning( "Invalid objtype for beta" ); value--; } PLA_Obj_global_length( beta, &size ); if ( size != 1 ){ PLA_Warning( "Invalid global length for beta" ); value--; } PLA_Obj_global_width( beta, &size ); if ( size != 1 ){ PLA_Warning( "Invalid global width for beta" ); value--; } /* Check if y is valid vector */ if ( y == NULL || !PLA_Valid_object( y ) ) { PLA_Warning( "Invalid object y" ); value--; } PLA_Obj_objtype( y, &objtype ); if ( objtype != PLA_MATRIX && objtype != PLA_MVECTOR && objtype != PLA_PMVECTOR ){ PLA_Warning( "Invalid objtype for y" ); value--; } PLA_Obj_project_onto( y, &proj_onto ); if ( proj_onto == PLA_PROJ_ONTO_COL ){ PLA_Obj_global_length( y, &length_y ); PLA_Obj_global_width( y, &width_y ); } else { PLA_Obj_global_width( y, &length_y ); PLA_Obj_global_length( y, &width_y ); } if ( width_y != 1 ){ PLA_Warning( "y is not of width 1" ); value--; } /* Check if dimensions match */ if ( transa == PLA_NO_TRANS || transa == PLA_CONJ ){ if ( length_A != length_y ){ PLA_Warning( "length of A does not match length of y" ); value--; } if ( width_A != length_x ){ PLA_Warning( "width of A does not match length of x" ); value--; } } else { if ( length_A != length_x ){ PLA_Warning( "length of A does not match length of x" ); value--; } if ( width_A != length_y ){ PLA_Warning( "width of A does not match length of y" ); value--; } } } if ( PLA_CHECK_AGAINST_SEQUENTIAL ){ PLA_Mscalar_create_conf_to( alpha, PLA_ALL_ROWS, PLA_ALL_COLS, &alpha_cpy ); PLA_Mscalar_create_conf_to( A, PLA_ALL_ROWS, PLA_ALL_COLS, &A_cpy ); PLA_Mscalar_create_conf_to( x, PLA_ALL_ROWS, PLA_ALL_COLS, &x_cpy ); PLA_Mscalar_create_conf_to( beta, PLA_ALL_ROWS, PLA_ALL_COLS, &beta_cpy ); PLA_Mscalar_create_conf_to( y, PLA_ALL_ROWS, PLA_ALL_COLS, &y_cpy ); PLA_Copy( alpha, alpha_cpy ); PLA_Copy( A, A_cpy ); PLA_Copy( x, x_cpy ); PLA_Copy( beta, beta_cpy ); if ( PLA_Local_equal_zero( beta_cpy ) ) PLA_Obj_set_to_zero( y_cpy ); else PLA_Copy( y, y_cpy ); } PLA_Routine_stack_pop( routine_name ); return value; } #define max(x,y) ( (x) > (y) ? (x) : (y) ) int PLA_Gemv_exit( int transa, PLA_Obj alpha, PLA_Obj A, PLA_Obj x, PLA_Obj beta, PLA_Obj y ) { int value = PLA_SUCCESS, size_malloced; double max_A, max_x, max_y, max_all, diff, PLA_Local_abs_max(), PLA_Local_abs_diff(); PLA_Obj one = NULL; char routine_name[ 35 ]; PLA_Routine_stack_push( "PLA_Gemv_exit" ); if ( PLA_CHECK_AGAINST_SEQUENTIAL ){ PLA_Obj y_temp = NULL; PLA_Mscalar_create_conf_to( y, PLA_ALL_ROWS, PLA_ALL_COLS, &y_temp ); PLA_Copy( y, y_temp ); max_A = PLA_Local_abs_max( A_cpy ); max_x = PLA_Local_abs_max( x_cpy ); PLA_Local_scal( beta_cpy, y_cpy ); max_y = PLA_Local_abs_max( y_cpy ); max_all = max( max_A, max( max_x, max_y ) ); PLA_Create_constants_conf_to( A, NULL, NULL, &one ); PLA_Local_gemv( transa, alpha_cpy, A_cpy, x_cpy, one, y_cpy ); diff = PLA_Local_abs_diff( y_temp, y_cpy ); if ( diff > 0.000001 * max_all ){ PLA_Warning( "PLA_Gemv: large relative error encountered" ); value--; } PLA_Obj_free( &alpha_cpy ); PLA_Obj_free( &A_cpy ); PLA_Obj_free( &x_cpy ); PLA_Obj_free( &beta_cpy ); PLA_Obj_free( &y_cpy ); PLA_Obj_free( &y_temp ); PLA_Obj_free( &one ); } size_malloced = PLA_Total_size_malloced( ); if ( size_malloced != old_size_malloced ) PLA_Warning( "PLA_Gemv: memory discrepency" ); PLA_Routine_stack_pop( routine_name ); if ( strcmp( routine_name, "PLA_Gemv_exit" ) != 0 ) PLA_Warning( "PLA_Gemv_exit: history stack corrupted" ); PLA_Routine_stack_pop( routine_name ); if ( strcmp( routine_name, "PLA_Gemv" ) != 0 ) PLA_Warning( "PLA_Gemv: history stack corrupted" ); return value; }