# include # include # include # include # include # include "padua.h" /******************************************************************************/ void filename_inc ( char *filename ) /******************************************************************************/ /* Purpose: FILENAME_INC increments a partially numeric file name. Discussion: It is assumed that the digits in the name, whether scattered or connected, represent a number that is to be increased by 1 on each call. If this number is all 9's on input, the output number is all 0's. Non-numeric letters of the name are unaffected. If the name is empty, then the routine stops. If the name contains no digits, the empty string is returned. Example: Input Output ----- ------ "a7to11.txt" "a7to12.txt" (typical case. Last digit incremented) "a7to99.txt" "a8to00.txt" (last digit incremented, with carry.) "a9to99.txt" "a0to00.txt" (wrap around) "cat.txt" " " (no digits to increment) " " STOP! (error) Licensing: This code is distributed under the MIT license. Modified: 22 November 2011 Author: John Burkardt Parameters: Input/output, char *FILENAME, the filename to be incremented. */ { int change; int n; char *t; n = strlen ( filename ); if ( n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FILENAME_INC - Fatal error!\n" ); fprintf ( stderr, " The input string is empty.\n" ); exit ( 1 ); } change = 0; t = filename + n - 1; while ( 0 < n ) { if ( '0' <= *t && *t <= '9' ) { change = change + 1; if ( *t == '9' ) { *t = '0'; } else { *t = *t + 1; return; } } t--; n--; } /* No digits were found. Return blank. */ if ( change == 0 ) { n = strlen ( filename ); t = filename + n - 1; while ( 0 < n ) { *t = ' '; t--; n--; } } return; } /******************************************************************************/ int padua_order ( int l ) /******************************************************************************/ /* Purpose: PADUA_ORDER returns the size of the Padua set of given level. Discussion: The Padua sets are indexed by a level that starts at 0. This function returns the number of points in each level. Example: Level Size ----- ---- 0 1 1 3 2 6 3 10 4 15 5 21 Licensing: This code is distributed under the MIT license. Modified: 09 February 2014 Author: John Burkardt Reference: Marco Caliari, Stefano de Marchi, Marco Vianello, Bivariate interpolation on the square at new nodal sets, Applied Mathematics and Computation, Volume 165, Number 2, 2005, pages 261-274. Parameters: Input, int L, the level of the set. 0 <= L Output, int PADUA_ORDER, the order (number of points) in the set. */ { int i; int n; n = 0; for ( i = 0; i <= l; i++ ) { n = n + ( l / 2 ) + 1; if ( ( l % 2 ) == 1 && ( i % 2 ) == 1 ) { n = n + 1; } } return n; } /******************************************************************************/ void padua_plot ( int l, char *filename ) /******************************************************************************/ /* Purpose: PADUA_PLOT plots the Padua points of given level. Licensing: This code is distributed under the MIT license. Modified: 09 June 2014 Author: John Burkardt Parameters: Input, int L, the level of the set. 0 <= L Input, char *FILENAME, a common filename prefix for the files to be created. */ { char command_filename[255]; FILE *command_unit; char data_filename[255]; FILE *data_unit; int j; int n; char plot_filename[255]; double *xy; n = padua_order ( l ); xy = padua_points ( l ); /* Create graphics data file. */ strcpy ( data_filename, filename ); strcat ( data_filename, "_data.txt" ); data_unit = fopen ( data_filename, "wt" ); for ( j = 0; j < n; j++ ) { fprintf ( data_unit, " %24.16g %24.16g\n", xy[0+j*2], xy[1+j*2] ); } fclose ( data_unit ); printf ( "\n" ); printf ( " Created data file '%s'.\n", data_filename ); /* Create graphics command file. */ strcpy ( command_filename, filename ); strcat ( command_filename, "_commands.txt" ); command_unit = fopen ( command_filename, "wt" ); fprintf ( command_unit, "# %s\n", command_filename ); fprintf ( command_unit, "#\n" ); fprintf ( command_unit, "# Usage:\n" ); fprintf ( command_unit, "# gnuplot < %s\n", command_filename ); fprintf ( command_unit, "#\n" ); fprintf ( command_unit, "set term png\n" ); strcpy ( plot_filename, filename ); strcat ( plot_filename, ".png" ); fprintf ( command_unit, "set output '%s'\n", plot_filename ); fprintf ( command_unit, "set xlabel '<--- X --->'\n" ); fprintf ( command_unit, "set ylabel '<--- Y --->'\n" ); fprintf ( command_unit, "set title 'Padua Points, Level %d\n", l ); fprintf ( command_unit, "set grid\n" ); fprintf ( command_unit, "set key off\n" ); fprintf ( command_unit, "set size ratio -1\n" ); fprintf ( command_unit, "set style data lines\n" ); fprintf ( command_unit, "set timestamp\n" ); fprintf ( command_unit, "plot [-1:+1] [-1:+1] '%s' using 1:2 with points lt 3 pt 3\n", data_filename ); fclose ( command_unit ); printf ( " Created command file '%s'.\n", command_filename ); /* Free memory. */ free ( xy ); return; } /******************************************************************************/ double *padua_points ( int l ) /******************************************************************************/ /* Purpose: PADUA_POINTS returns the Padua points of order N. Licensing: This code is distributed under the MIT license. Modified: 09 June 2014 Author: John Burkardt Reference: Marco Caliari, Stefano de Marchi, Marco Vianello, Bivariate interpolation on the square at new nodal sets, Applied Mathematics and Computation, Volume 165, Number 2, 2005, pages 261-274. Parameters: Input, int L, the level of the set. 0 <= L Output, double PADUA_POINTS[2*(((L+1)*(L+2))/2)], the Padua points. */ { double angle1; double angle2; int i; int j; int j_hi; int k; int n; const double r8_pi = 3.141592653589793; double *xy; n = ( ( l + 1 ) * ( l + 2 ) ) / 2; xy = ( double * ) malloc ( 2 * n * sizeof ( double ) ); if ( l == 0 ) { xy[0+0*2] = 0.0; xy[1+0*2] = 0.0; return xy; } k = 0; for ( i = 0; i <= l; i++ ) { j_hi = ( l / 2 ) + 1; if ( ( l % 2 ) == 1 && ( i % 2 ) == 1 ) { j_hi = j_hi + 1; } for ( j = 1; j <= j_hi; j++ ) { if ( i * 2 == l ) { xy[0+k*2] = 0.0; } else { angle1 = ( double ) ( i ) * r8_pi / ( double ) ( l ); xy[0+k*2] = cos ( angle1 ); } if ( ( i % 2 ) == 0 ) { if ( 2 * ( 2 * j - 1 ) == l + 1 ) { xy[1+k*2] = 0.0; } else { angle2 = ( double ) ( 2 * j - 1 ) * r8_pi / ( double ) ( l + 1 ); xy[1+k*2] = cos ( angle2 ); } } else { if ( 2 * ( 2 * j - 2 ) == l + 1 ) { xy[1+k*2] = 0.0; } else { angle2 = ( double ) ( 2 * j - 2 ) * r8_pi / ( double ) ( l + 1 ); xy[1+k*2] = cos ( angle2 ); } } k = k + 1; } } return xy; } /******************************************************************************/ double *padua_points_set ( int l ) /******************************************************************************/ /* Purpose: PADUA_POINTS_SET sets the Padua points. Licensing: This code is distributed under the MIT license. Modified: 30 August 2016 Author: John Burkardt Reference: Marco Caliari, Stefano de Marchi, Marco Vianello, Bivariate interpolation on the square at new nodal sets, Applied Mathematics and Computation, Volume 165, Number 2, 2005, pages 261-274. Parameters: Input, int L, the level. 0 <= L <= 10. Output, double PADUA_POINTS_SET[2*N], the Padua points. */ { int j1; int j1_hi; int j2; int n; double t; double *xy; double xy00[2*1] = { 0.000000000000000, 0.000000000000000 }; double xy01[2*3] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, 1.000000000000000, 1.000000000000000, 0.000000000000000 }; double xy02[2*6] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, 0.5000000000000001, 0.000000000000000, -0.4999999999999998, 0.000000000000000, 1.000000000000000, 1.000000000000000, -1.000000000000000, 1.000000000000000, 0.5000000000000001 }; double xy03[2*10] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, 0.000000000000000, -1.000000000000000, 1.000000000000000, -0.4999999999999998, -0.7071067811865475, -0.4999999999999998, 0.7071067811865476, 0.5000000000000001, -1.000000000000000, 0.5000000000000001, 0.000000000000000, 0.5000000000000001, 1.000000000000000, 1.000000000000000, -0.7071067811865475, 1.000000000000000, 0.7071067811865476 }; double xy04[2*15] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, -0.3090169943749473, -1.000000000000000, 0.8090169943749475, -0.7071067811865475, -0.8090169943749473, -0.7071067811865475, 0.3090169943749475, -0.7071067811865475, 1.000000000000000, 0.000000000000000, -1.000000000000000, 0.000000000000000, -0.3090169943749473, 0.000000000000000, 0.8090169943749475, 0.7071067811865476, -0.8090169943749473, 0.7071067811865476, 0.3090169943749475, 0.7071067811865476, 1.000000000000000, 1.000000000000000, -1.000000000000000, 1.000000000000000, -0.3090169943749473, 1.000000000000000, 0.8090169943749475 }; double xy05[2*21] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, -0.4999999999999998, -1.000000000000000, 0.5000000000000001, -1.000000000000000, 1.000000000000000, -0.8090169943749473, -0.8660254037844387, -0.8090169943749473, 0.000000000000000, -0.8090169943749473, 0.8660254037844387, -0.3090169943749473, -1.000000000000000, -0.3090169943749473, -0.4999999999999998, -0.3090169943749473, 0.5000000000000001, -0.3090169943749473, 1.000000000000000, 0.3090169943749475, -0.8660254037844387, 0.3090169943749475, 0.000000000000000, 0.3090169943749475, 0.8660254037844387, 0.8090169943749475, -1.000000000000000, 0.8090169943749475, -0.4999999999999998, 0.8090169943749475, 0.5000000000000001, 0.8090169943749475, 1.000000000000000, 1.000000000000000, -0.8660254037844387, 1.000000000000000, 0.000000000000000, 1.000000000000000, 0.8660254037844387 }; double xy06[2*28] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, -0.6234898018587335, -1.000000000000000, 0.2225209339563144, -1.000000000000000, 0.9009688679024191, -0.8660254037844387, -0.9009688679024190, -0.8660254037844387, -0.2225209339563143, -0.8660254037844387, 0.6234898018587336, -0.8660254037844387, 1.000000000000000, -0.4999999999999998, -1.000000000000000, -0.4999999999999998, -0.6234898018587335, -0.4999999999999998, 0.2225209339563144, -0.4999999999999998, 0.9009688679024191, 0.000000000000000, -0.9009688679024190, 0.000000000000000, -0.2225209339563143, 0.000000000000000, 0.6234898018587336, 0.000000000000000, 1.000000000000000, 0.5000000000000001, -1.000000000000000, 0.5000000000000001, -0.6234898018587335, 0.5000000000000001, 0.2225209339563144, 0.5000000000000001, 0.9009688679024191, 0.8660254037844387, -0.9009688679024190, 0.8660254037844387, -0.2225209339563143, 0.8660254037844387, 0.6234898018587336, 0.8660254037844387, 1.000000000000000, 1.000000000000000, -1.000000000000000, 1.000000000000000, -0.6234898018587335, 1.000000000000000, 0.2225209339563144, 1.000000000000000, 0.9009688679024191 }; double xy07[2*36] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, -0.7071067811865475, -1.000000000000000, 0.000000000000000, -1.000000000000000, 0.7071067811865476, -1.000000000000000, 1.000000000000000, -0.9009688679024190, -0.9238795325112867, -0.9009688679024190, -0.3826834323650897, -0.9009688679024190, 0.3826834323650898, -0.9009688679024190, 0.9238795325112867, -0.6234898018587335, -1.000000000000000, -0.6234898018587335, -0.7071067811865475, -0.6234898018587335, 0.000000000000000, -0.6234898018587335, 0.7071067811865476, -0.6234898018587335, 1.000000000000000, -0.2225209339563143, -0.9238795325112867, -0.2225209339563143, -0.3826834323650897, -0.2225209339563143, 0.3826834323650898, -0.2225209339563143, 0.9238795325112867, 0.2225209339563144, -1.000000000000000, 0.2225209339563144, -0.7071067811865475, 0.2225209339563144, 0.000000000000000, 0.2225209339563144, 0.7071067811865476, 0.2225209339563144, 1.000000000000000, 0.6234898018587336, -0.9238795325112867, 0.6234898018587336, -0.3826834323650897, 0.6234898018587336, 0.3826834323650898, 0.6234898018587336, 0.9238795325112867, 0.9009688679024191, -1.000000000000000, 0.9009688679024191, -0.7071067811865475, 0.9009688679024191, 0.000000000000000, 0.9009688679024191, 0.7071067811865476, 0.9009688679024191, 1.000000000000000, 1.000000000000000, -0.9238795325112867, 1.000000000000000, -0.3826834323650897, 1.000000000000000, 0.3826834323650898, 1.000000000000000, 0.9238795325112867 }; double xy08[2*45] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, -0.7660444431189779, -1.000000000000000, -0.1736481776669303, -1.000000000000000, 0.5000000000000001, -1.000000000000000, 0.9396926207859084, -0.9238795325112867, -0.9396926207859083, -0.9238795325112867, -0.4999999999999998, -0.9238795325112867, 0.1736481776669304, -0.9238795325112867, 0.7660444431189780, -0.9238795325112867, 1.000000000000000, -0.7071067811865475, -1.000000000000000, -0.7071067811865475, -0.7660444431189779, -0.7071067811865475, -0.1736481776669303, -0.7071067811865475, 0.5000000000000001, -0.7071067811865475, 0.9396926207859084, -0.3826834323650897, -0.9396926207859083, -0.3826834323650897, -0.4999999999999998, -0.3826834323650897, 0.1736481776669304, -0.3826834323650897, 0.7660444431189780, -0.3826834323650897, 1.000000000000000, 0.000000000000000, -1.000000000000000, 0.000000000000000, -0.7660444431189779, 0.000000000000000, -0.1736481776669303, 0.000000000000000, 0.5000000000000001, 0.000000000000000, 0.9396926207859084, 0.3826834323650898, -0.9396926207859083, 0.3826834323650898, -0.4999999999999998, 0.3826834323650898, 0.1736481776669304, 0.3826834323650898, 0.7660444431189780, 0.3826834323650898, 1.000000000000000, 0.7071067811865476, -1.000000000000000, 0.7071067811865476, -0.7660444431189779, 0.7071067811865476, -0.1736481776669303, 0.7071067811865476, 0.5000000000000001, 0.7071067811865476, 0.9396926207859084, 0.9238795325112867, -0.9396926207859083, 0.9238795325112867, -0.4999999999999998, 0.9238795325112867, 0.1736481776669304, 0.9238795325112867, 0.7660444431189780, 0.9238795325112867, 1.000000000000000, 1.000000000000000, -1.000000000000000, 1.000000000000000, -0.7660444431189779, 1.000000000000000, -0.1736481776669303, 1.000000000000000, 0.5000000000000001, 1.000000000000000, 0.9396926207859084 }; double xy09[2*55] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, -0.8090169943749473, -1.000000000000000, -0.3090169943749473, -1.000000000000000, 0.3090169943749475, -1.000000000000000, 0.8090169943749475, -1.000000000000000, 1.000000000000000, -0.9396926207859083, -0.9510565162951535, -0.9396926207859083, -0.5877852522924730, -0.9396926207859083, 0.000000000000000, -0.9396926207859083, 0.5877852522924731, -0.9396926207859083, 0.9510565162951535, -0.7660444431189779, -1.000000000000000, -0.7660444431189779, -0.8090169943749473, -0.7660444431189779, -0.3090169943749473, -0.7660444431189779, 0.3090169943749475, -0.7660444431189779, 0.8090169943749475, -0.7660444431189779, 1.000000000000000, -0.4999999999999998, -0.9510565162951535, -0.4999999999999998, -0.5877852522924730, -0.4999999999999998, 0.000000000000000, -0.4999999999999998, 0.5877852522924731, -0.4999999999999998, 0.9510565162951535, -0.1736481776669303, -1.000000000000000, -0.1736481776669303, -0.8090169943749473, -0.1736481776669303, -0.3090169943749473, -0.1736481776669303, 0.3090169943749475, -0.1736481776669303, 0.8090169943749475, -0.1736481776669303, 1.000000000000000, 0.1736481776669304, -0.9510565162951535, 0.1736481776669304, -0.5877852522924730, 0.1736481776669304, 0.000000000000000, 0.1736481776669304, 0.5877852522924731, 0.1736481776669304, 0.9510565162951535, 0.5000000000000001, -1.000000000000000, 0.5000000000000001, -0.8090169943749473, 0.5000000000000001, -0.3090169943749473, 0.5000000000000001, 0.3090169943749475, 0.5000000000000001, 0.8090169943749475, 0.5000000000000001, 1.000000000000000, 0.7660444431189780, -0.9510565162951535, 0.7660444431189780, -0.5877852522924730, 0.7660444431189780, 0.000000000000000, 0.7660444431189780, 0.5877852522924731, 0.7660444431189780, 0.9510565162951535, 0.9396926207859084, -1.000000000000000, 0.9396926207859084, -0.8090169943749473, 0.9396926207859084, -0.3090169943749473, 0.9396926207859084, 0.3090169943749475, 0.9396926207859084, 0.8090169943749475, 0.9396926207859084, 1.000000000000000, 1.000000000000000, -0.9510565162951535, 1.000000000000000, -0.5877852522924730, 1.000000000000000, 0.000000000000000, 1.000000000000000, 0.5877852522924731, 1.000000000000000, 0.9510565162951535 }; double xy10[2*66] = { -1.000000000000000, -1.000000000000000, -1.000000000000000, -0.8412535328311811, -1.000000000000000, -0.4154150130018863, -1.000000000000000, 0.1423148382732851, -1.000000000000000, 0.6548607339452851, -1.000000000000000, 0.9594929736144974, -0.9510565162951535, -0.9594929736144974, -0.9510565162951535, -0.6548607339452850, -0.9510565162951535, -0.1423148382732850, -0.9510565162951535, 0.4154150130018864, -0.9510565162951535, 0.8412535328311812, -0.9510565162951535, 1.000000000000000, -0.8090169943749473, -1.000000000000000, -0.8090169943749473, -0.8412535328311811, -0.8090169943749473, -0.4154150130018863, -0.8090169943749473, 0.1423148382732851, -0.8090169943749473, 0.6548607339452851, -0.8090169943749473, 0.9594929736144974, -0.5877852522924730, -0.9594929736144974, -0.5877852522924730, -0.6548607339452850, -0.5877852522924730, -0.1423148382732850, -0.5877852522924730, 0.4154150130018864, -0.5877852522924730, 0.8412535328311812, -0.5877852522924730, 1.000000000000000, -0.3090169943749473, -1.000000000000000, -0.3090169943749473, -0.8412535328311811, -0.3090169943749473, -0.4154150130018863, -0.3090169943749473, 0.1423148382732851, -0.3090169943749473, 0.6548607339452851, -0.3090169943749473, 0.9594929736144974, 0.000000000000000, -0.9594929736144974, 0.000000000000000, -0.6548607339452850, 0.000000000000000, -0.1423148382732850, 0.000000000000000, 0.4154150130018864, 0.000000000000000, 0.8412535328311812, 0.000000000000000, 1.000000000000000, 0.3090169943749475, -1.000000000000000, 0.3090169943749475, -0.8412535328311811, 0.3090169943749475, -0.4154150130018863, 0.3090169943749475, 0.1423148382732851, 0.3090169943749475, 0.6548607339452851, 0.3090169943749475, 0.9594929736144974, 0.5877852522924731, -0.9594929736144974, 0.5877852522924731, -0.6548607339452850, 0.5877852522924731, -0.1423148382732850, 0.5877852522924731, 0.4154150130018864, 0.5877852522924731, 0.8412535328311812, 0.5877852522924731, 1.000000000000000, 0.8090169943749475, -1.000000000000000, 0.8090169943749475, -0.8412535328311811, 0.8090169943749475, -0.4154150130018863, 0.8090169943749475, 0.1423148382732851, 0.8090169943749475, 0.6548607339452851, 0.8090169943749475, 0.9594929736144974, 0.9510565162951535, -0.9594929736144974, 0.9510565162951535, -0.6548607339452850, 0.9510565162951535, -0.1423148382732850, 0.9510565162951535, 0.4154150130018864, 0.9510565162951535, 0.8412535328311812, 0.9510565162951535, 1.000000000000000, 1.000000000000000, -1.000000000000000, 1.000000000000000, -0.8412535328311811, 1.000000000000000, -0.4154150130018863, 1.000000000000000, 0.1423148382732851, 1.000000000000000, 0.6548607339452851, 1.000000000000000, 0.9594929736144974 }; n = ( ( l + 1 ) * ( l + 2 ) ) / 2; if ( l == 0 ) { xy = r8vec_copy_new ( 2 * n, xy00 ); } else if ( l == 1 ) { xy = r8vec_copy_new ( 2 * n, xy01 ); } else if ( l == 2 ) { xy = r8vec_copy_new ( 2 * n, xy02 ); } else if ( l == 3 ) { xy = r8vec_copy_new ( 2 * n, xy03 ); } else if ( l == 4 ) { xy = r8vec_copy_new ( 2 * n, xy04 ); } else if ( l == 5 ) { xy = r8vec_copy_new ( 2 * n, xy05 ); } else if ( l == 6 ) { xy = r8vec_copy_new ( 2 * n, xy06 ); } else if ( l == 7 ) { xy = r8vec_copy_new ( 2 * n, xy07 ); } else if ( l == 8 ) { xy = r8vec_copy_new ( 2 * n, xy08 ); } else if ( l == 9 ) { xy = r8vec_copy_new ( 2 * n, xy09 ); } else if ( l == 10 ) { xy = r8vec_copy_new ( 2 * n, xy10 ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "PADUA_POINTS_SET - Fatal error!\n" ); fprintf ( stderr, " Illegal value of L = %d\n", l ); fprintf ( stderr, " Legal values are 1 through 10.\n" ); exit ( 1 ); } /* Reverse data to match published information. */ j1_hi = ( n - 1 ) / 2; for ( j1 = 0; j1 < j1_hi; j1++ ) { j2 = n - 1 - j1; t = xy[0+2*j1]; xy[0+2*j1] = xy[0+2*j2]; xy[0+2*j2] = t; t = xy[1+2*j1]; xy[1+2*j1] = xy[1+2*j2]; xy[1+2*j2] = t; } return xy; } /******************************************************************************/ double *padua_weights_set ( int l ) /******************************************************************************/ /* Purpose: PADUA_WEIGHTS_SET sets quadrature weights for the Padua points. Licensing: This code is distributed under the MIT license. Modified: 11 June 2014 Author: John Burkardt Reference: Marco Caliari, Stefano de Marchi, Marco Vianello, Bivariate interpolation on the square at new nodal sets, Applied Mathematics and Computation, Volume 165, Number 2, 2005, pages 261-274. Parameters: Input, int L, the level. 0 <= L <= 10. Output, double W[((L+1)*(L+2))/2], the quadrature weights. */ { int n; double *w; n = ( ( l + 1 ) * ( l + 2 ) ) / 2; w = ( double * ) malloc ( n * sizeof ( double ) ); if ( l == 0 ) { w[ 0] = 4.000000000000000E+00; } else if ( l == 1 ) { w[ 0] = 1.000000000000000E+00; w[ 1] = 1.000000000000000E+00; w[ 2] = 2.000000000000000E+00; } else if ( l == 2 ) { w[ 0] = 0.0E+00; w[ 1] = 0.6666666666666663E+00; w[ 2] = 2.222222222222222E+00; w[ 3] = 0.4444444444444444E+00; w[ 4] = 0.0E+00; w[ 5] = 0.6666666666666664E+00; } else if ( l == 3 ) { w[ 0] = -0.5555555555555480E-01; w[ 1] = 0.3333333333333331E+00; w[ 2] = -0.5555555555555580E-01; w[ 3] = 0.8888888888888886E+00; w[ 4] = 0.8888888888888893E+00; w[ 5] = 0.2222222222222224E+00; w[ 6] = 1.333333333333333E+00; w[ 7] = 0.2222222222222220E+00; w[ 8] = 0.1111111111111109E+00; w[ 9] = 0.1111111111111112E+00; } else if ( l == 4 ) { w[ 0] = -0.8888888888888932E-02; w[ 1] = 0.8104919101110961E-01; w[ 2] = 0.6117303121111219E-01; w[ 3] = 0.3874097078666789E+00; w[ 4] = 0.6259236254666545E+00; w[ 5] = 0.5333333333333362E-01; w[ 6] = 0.7111111111111067E-01; w[ 7] = 0.9830822022444241E+00; w[ 8] = 0.5458066866444642E+00; w[ 9] = 0.3874097078666780E+00; w[10] = 0.6259236254666568E+00; w[11] = 0.5333333333333383E-01; w[12] = -0.8888888888888703E-02; w[13] = 0.8104919101110968E-01; w[14] = 0.6117303121111135E-01; } else if ( l == 5 ) { w[ 0] = -0.1037037037037093E-01; w[ 1] = 0.5037037037036911E-01; w[ 2] = 0.5037037037037081E-01; w[ 3] = -0.1037037037036947E-01; w[ 4] = 0.1876963678740801E+00; w[ 5] = 0.3460933466518654E+00; w[ 6] = 0.1876963678740763E+00; w[ 7] = 0.4514390511851724E-01; w[ 8] = 0.5541130536814713E+00; w[ 9] = 0.5541130536814728E+00; w[10] = 0.4514390511851834E-01; w[11] = 0.2804517802740705E+00; w[12] = 0.6376103570518378E+00; w[13] = 0.2804517802740683E+00; w[14] = 0.3189313191851883E-01; w[15] = 0.3288499092814910E+00; w[16] = 0.3288499092814925E+00; w[17] = 0.3189313191851956E-01; w[18] = 0.2074074074074123E-01; w[19] = 0.3851851851851849E-01; w[20] = 0.2074074074074051E-01; } else if ( l == 6 ) { w[ 0] = -0.3023431594858565E-02; w[ 1] = 0.1957267632451884E-01; w[ 2] = 0.2633929313290840E-01; w[ 3] = 0.1425431928029237E-01; w[ 4] = 0.1006383046329639E+00; w[ 5] = 0.2208900184526934E+00; w[ 6] = 0.1743144584714012E+00; w[ 7] = 0.1209372637943976E-01; w[ 8] = 0.1934996220710680E-01; w[ 9] = 0.3245064820875231E+00; w[10] = 0.4027058473592984E+00; w[11] = 0.1677234226317961E+00; w[12] = 0.1953319357827178E+00; w[13] = 0.4489633053035124E+00; w[14] = 0.3721824611057551E+00; w[15] = 0.2479213907785274E-01; w[16] = 0.1934996220710561E-01; w[17] = 0.3245064820875153E+00; w[18] = 0.4027058473592959E+00; w[19] = 0.1677234226317944E+00; w[20] = 0.1006383046329745E+00; w[21] = 0.2208900184526933E+00; w[22] = 0.1743144584714027E+00; w[23] = 0.1209372637944051E-01; w[24] = -0.3023431594861990E-02; w[25] = 0.1957267632451757E-01; w[26] = 0.2633929313290797E-01; w[27] = 0.1425431928029198E-01; } else if ( l == 7 ) { w[ 0] = -0.3287981859413765E-02; w[ 1] = 0.1337868480725671E-01; w[ 2] = 0.2063492063491996E-01; w[ 3] = 0.1337868480725546E-01; w[ 4] = -0.3287981859408898E-02; w[ 5] = 0.5949324721885513E-01; w[ 6] = 0.1306477599993571E+00; w[ 7] = 0.1306477599993581E+00; w[ 8] = 0.5949324721885061E-01; w[ 9] = 0.1263869091685831E-01; w[10] = 0.1979944935601103E+00; w[11] = 0.2832184784823740E+00; w[12] = 0.1979944935601143E+00; w[13] = 0.1263869091685747E-01; w[14] = 0.1221817987389771E+00; w[15] = 0.3150266070593529E+00; w[16] = 0.3150266070593440E+00; w[17] = 0.1221817987389802E+00; w[18] = 0.1771365352315134E-01; w[19] = 0.2490926964598258E+00; w[20] = 0.3408041116306980E+00; w[21] = 0.2490926964598291E+00; w[22] = 0.1771365352314976E-01; w[23] = 0.9646986307476696E-01; w[24] = 0.2557725606433917E+00; w[25] = 0.2557725606433927E+00; w[26] = 0.9646986307476431E-01; w[27] = 0.8649923133686802E-02; w[28] = 0.1062007918394705E+00; w[29] = 0.1505805844901012E+00; w[30] = 0.1062007918394705E+00; w[31] = 0.8649923133690016E-02; w[32] = 0.6355881462931014E-02; w[33] = 0.1405228180237514E-01; w[34] = 0.1405228180237651E-01; w[35] = 0.6355881462928496E-02; } else if ( l == 8 ) { w[ 0] = -0.1269841269835311E-02; w[ 1] = 0.6706089639041270E-02; w[ 2] = 0.1111455441352989E-01; w[ 3] = 0.1026455026455282E-01; w[ 4] = 0.4930678698742625E-02; w[ 5] = 0.3633146869162523E-01; w[ 6] = 0.8838322767333079E-01; w[ 7] = 0.9965911758463214E-01; w[ 8] = 0.6400185533755555E-01; w[ 9] = 0.4061629144893127E-02; w[10] = 0.6772486772485166E-02; w[11] = 0.1258344472781388E+00; w[12] = 0.1927501398511116E+00; w[13] = 0.1699470899470907E+00; w[14] = 0.6342599488133535E-01; w[15] = 0.8376332474107638E-01; w[16] = 0.2170841444607031E+00; w[17] = 0.2477307250801775E+00; w[18] = 0.1648098048612226E+00; w[19] = 0.1004771829779292E-01; w[20] = 0.1015873015872910E-01; w[21] = 0.1784328991205164E+00; w[22] = 0.2729409493576765E+00; w[23] = 0.2364021164021134E+00; w[24] = 0.8936689226256009E-01; w[25] = 0.8376332474107701E-01; w[26] = 0.2170841444607054E+00; w[27] = 0.2477307250801761E+00; w[28] = 0.1648098048612200E+00; w[29] = 0.1004771829779330E-01; w[30] = 0.6772486772485237E-02; w[31] = 0.1258344472781358E+00; w[32] = 0.1927501398511135E+00; w[33] = 0.1699470899470926E+00; w[34] = 0.6342599488133838E-01; w[35] = 0.3633146869162453E-01; w[36] = 0.8838322767332588E-01; w[37] = 0.9965911758463601E-01; w[38] = 0.6400185533755502E-01; w[39] = 0.4061629144888279E-02; w[40] = -0.1269841269836355E-02; w[41] = 0.6706089639046927E-02; w[42] = 0.1111455441352761E-01; w[43] = 0.1026455026454956E-01; w[44] = 0.4930678698747173E-02; } else if ( l == 9 ) { w[ 0] = -0.1368606701945113E-02; w[ 1] = 0.4837977417140975E-02; w[ 2] = 0.8876308297144902E-02; w[ 3] = 0.8876308297143068E-02; w[ 4] = 0.4837977417150492E-02; w[ 5] = -0.1368606701935084E-02; w[ 6] = 0.2425285860992349E-01; w[ 7] = 0.5727330842923516E-01; w[ 8] = 0.7008257906578071E-01; w[ 9] = 0.5727330842922034E-01; w[10] = 0.2425285860989794E-01; w[11] = 0.4659404339099723E-02; w[12] = 0.8354521980498550E-01; w[13] = 0.1370796991940044E+00; w[14] = 0.1370796991940248E+00; w[15] = 0.8354521980500107E-01; w[16] = 0.4659404339109654E-02; w[17] = 0.5564545640233619E-01; w[18] = 0.1524391996823315E+00; w[19] = 0.1877107583774149E+00; w[20] = 0.1524391996823176E+00; w[21] = 0.5564545640232402E-01; w[22] = 0.8186176158691754E-02; w[23] = 0.1295355639606716E+00; w[24] = 0.2061407656847711E+00; w[25] = 0.2061407656847630E+00; w[26] = 0.1295355639606894E+00; w[27] = 0.8186176158692687E-02; w[28] = 0.6234969028097752E-01; w[29] = 0.1730419031522391E+00; w[30] = 0.2169418247419051E+00; w[31] = 0.1730419031522361E+00; w[32] = 0.6234969028097048E-01; w[33] = 0.7506172839505762E-02; w[34] = 0.1142161960569350E+00; w[35] = 0.1802176663769002E+00; w[36] = 0.1802176663769038E+00; w[37] = 0.1142161960569279E+00; w[38] = 0.7506172839512260E-02; w[39] = 0.4031900987631698E-01; w[40] = 0.1142976211857364E+00; w[41] = 0.1413353845521477E+00; w[42] = 0.1142976211857414E+00; w[43] = 0.4031900987631700E-01; w[44] = 0.3239075586856897E-02; w[45] = 0.4317587564913915E-01; w[46] = 0.7015250533601934E-01; w[47] = 0.7015250533601930E-01; w[48] = 0.4317587564913908E-01; w[49] = 0.3239075586852207E-02; w[50] = 0.2550690557469151E-02; w[51] = 0.6084230077461027E-02; w[52] = 0.7421516754852508E-02; w[53] = 0.6084230077458821E-02; w[54] = 0.2550690557473353E-02; } else if ( l == 10 ) { w[ 0] = -0.6240762604463766E-03; w[ 1] = 0.2843227149025789E-02; w[ 2] = 0.5250031948150784E-02; w[ 3] = 0.5891746241568810E-02; w[ 4] = 0.4705736485964679E-02; w[ 5] = 0.2135354637732944E-02; w[ 6] = 0.1610939653924566E-01; w[ 7] = 0.4099595211758227E-01; w[ 8] = 0.5326500934654063E-01; w[ 9] = 0.4863338516658277E-01; w[10] = 0.2843474741781434E-01; w[11] = 0.1719619179693151E-02; w[12] = 0.2883769745121509E-02; w[13] = 0.5724711668876453E-01; w[14] = 0.9659872841640438E-01; w[15] = 0.1053210323353631E+00; w[16] = 0.8066212502628711E-01; w[17] = 0.2855765663647366E-01; w[18] = 0.3981286043310814E-01; w[19] = 0.1090390674981577E+00; w[20] = 0.1430169021081585E+00; w[21] = 0.1313686303763064E+00; w[22] = 0.7932850918298831E-01; w[23] = 0.4610696968783255E-02; w[24] = 0.5086495679684716E-02; w[25] = 0.9311356395361167E-01; w[26] = 0.1562320334111262E+00; w[27] = 0.1696057154254139E+00; w[28] = 0.1283581371975154E+00; w[29] = 0.4603059518094556E-01; w[30] = 0.4894888812994630E-01; w[31] = 0.1347281473526573E+00; w[32] = 0.1764193542601264E+00; w[33] = 0.1635037456303485E+00; w[34] = 0.9822749154565460E-01; w[35] = 0.5704840613923174E-02; w[36] = 0.5086495679679268E-02; w[37] = 0.9311356395362781E-01; w[38] = 0.1562320334111511E+00; w[39] = 0.1696057154253968E+00; w[40] = 0.1283581371975113E+00; w[41] = 0.4603059518094044E-01; w[42] = 0.3981286043311782E-01; w[43] = 0.1090390674981293E+00; w[44] = 0.1430169021081508E+00; w[45] = 0.1313686303763217E+00; w[46] = 0.7932850918299997E-01; w[47] = 0.4610696968790496E-02; w[48] = 0.2883769745110260E-02; w[49] = 0.5724711668875122E-01; w[50] = 0.9659872841642343E-01; w[51] = 0.1053210323353932E+00; w[52] = 0.8066212502626474E-01; w[53] = 0.2855765663644533E-01; w[54] = 0.1610939653928420E-01; w[55] = 0.4099595211758404E-01; w[56] = 0.5326500934649123E-01; w[57] = 0.4863338516656233E-01; w[58] = 0.2843474741784810E-01; w[59] = 0.1719619179720036E-02; w[60] = -0.6240762604606350E-03; w[61] = 0.2843227149011163E-02; w[62] = 0.5250031948172295E-02; w[63] = 0.5891746241587802E-02; w[64] = 0.4705736485965663E-02; w[65] = 0.2135354637703863E-02; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "PADUA_WEIGHTS_SET - Fatal error\n" ); fprintf ( stderr, " Illegal value of L = %d\n", l ); fprintf ( stderr, " Legal values are 0 through 10.\n" ); exit ( 1 ); } /* Reverse data to match published information. */ r8vec_reverse ( n, w ); return w; } /******************************************************************************/ double *padua_weights ( int l ) /******************************************************************************/ /* Purpose: PADUA_WEIGHTS returns quadrature weights do Padua points. Discussion: The order of the weights corresponds to the ordering used by the companion function padua_points(). Caliari, de Marchi and Vianello supplied a MATLAB code pdwtsMM which carries out this same computation in a way that makes more efficient use of MATLAB's vector and matrix capabilities. This version of the computation was painfully rewritten to display the individual scalar computations, so that it could be translated into other languages. Licensing: This code is distributed under the MIT license. Modified: 11 June 2014 Author: John Burkardt Reference: Marco Caliari, Stefano de Marchi, Marco Vianello, Bivariate interpolation on the square at new nodal sets, Applied Mathematics and Computation, Volume 165, Number 2, 2005, pages 261-274. Parameters: Input, int L, the level of the set. 0 <= L Output, double PADUA_WEIGHTS[(L+1)*(L+2)/2], the quadrature weights. */ { double angle; int i; int i2; int j; int j2; int lp1h; int lp2h; int lp3h; double mi; double mj; double *mom; int n; const double r8_pi = 3.141592653589793; double *te1; double *te2; double *tmteo; double *tmtoe; double *to1; double *to2; double *w; double *w1; double *w2; n = ( ( l + 1 ) * ( l + 2 ) ) / 2; w = ( double * ) malloc ( n * sizeof ( double ) ); if ( l == 0 ) { w[0] = 4.0; return w; } /* Relatives of L/2: */ lp1h = ( l + 1 ) / 2; lp2h = ( l + 2 ) / 2; lp3h = ( l + 3 ) / 2; /* TE1, TE2, TO1, TO2: Even and odd Chebyshev polynomials on subgrids 1 and 2. */ te1 = ( double * ) malloc ( lp2h * lp2h * sizeof ( double ) ); for ( j = 0; j < lp2h; j++ ) { for ( i = 0; i < lp2h; i++ ) { angle = r8_pi * ( double ) ( 2 * i * 2 * j ) / ( double ) ( l ); te1[i+j*lp2h] = cos ( angle ); } } for ( j = 0; j < lp2h; j++ ) { for ( i = 1; i < lp2h; i++ ) { te1[i+j*lp2h] = te1[i+j*lp2h] * sqrt ( 2.0 ); } } to1 = ( double * ) malloc ( lp2h * lp1h * sizeof ( double ) ); for ( j = 0; j < lp1h; j++ ) { for ( i = 0; i < lp2h; i++ ) { angle = r8_pi * ( double ) ( 2 * i * ( 2 * j + 1 ) ) / ( double ) ( l ); to1[i+j*lp2h] = cos ( angle ); } } for ( j = 0; j < lp1h; j++ ) { for ( i = 1; i < lp2h; i++ ) { to1[i+j*lp2h] = to1[i+j*lp2h] * sqrt ( 2.0 ); } } te2 = ( double * ) malloc ( lp2h * lp3h * sizeof ( double ) ); for ( j = 0; j < lp3h; j++ ) { for ( i = 0; i < lp2h; i++ ) { angle = r8_pi * ( double ) ( 2 * i * 2 * j ) / ( double ) ( l + 1 ); te2[i+j*lp2h] = cos ( angle ); } } for ( j = 0; j < lp3h; j++ ) { for ( i = 1; i < lp2h; i++ ) { te2[i+j*lp2h] = te2[i+j*lp2h] * sqrt ( 2.0 ); } } to2 = ( double * ) malloc ( lp2h * lp2h * sizeof ( double ) ); for ( j = 0; j < lp2h; j++ ) { for ( i = 0; i < lp2h; i++ ) { angle = r8_pi * ( double ) ( 2 * i * ( 2 * j + 1 ) ) / ( double ) ( l + 1 ); to2[i+j*lp2h] = cos ( angle ); } } for ( j = 0; j < lp2h; j++ ) { for ( i = 1; i < lp2h; i++ ) { to2[i+j*lp2h] = to2[i+j*lp2h] * sqrt ( 2.0 ); } } /* MOM: Moments matrix do even * even pairs. */ mom = ( double * ) malloc ( lp2h * lp2h * sizeof ( double ) ); for ( j = 0; j < lp2h; j++ ) { mj = 2.0 * sqrt ( 2.0 ) / ( double ) ( 1 - pow ( 2 * j, 2 ) ); for ( i = 0; i < lp2h - j; i++ ) { mi = 2.0 * sqrt ( 2.0 ) / ( double ) ( 1 - pow ( 2 * i, 2 ) ); mom[i+j*lp2h] = mi * mj; } } i = 0; for ( j = 0; j < lp2h; j++ ) { mom[i+j*lp2h] = mom[i+j*lp2h] / sqrt ( 2.0 ); } j = 0; for ( i = 0; i < lp2h; i++ ) { mom[i+j*lp2h] = mom[i+j*lp2h] / sqrt ( 2.0 ); } if ( ( l % 2 ) == 0 ) { i = lp2h-1; j = 0; mom[i+j*lp2h] = mom[i+j*lp2h] / 2.0; } /* TMTOE and TMTEO: matrix products. */ tmtoe = ( double * ) malloc ( lp2h * lp2h * sizeof ( double ) ); for ( j = 0; j < lp2h; j++ ) { for ( i = 0; i < lp2h; i++ ) { tmtoe[i+j*lp2h] = 0.0; } } for ( j2 = 0; j2 < lp2h; j2++ ) { for ( i2 = 0; i2 < lp2h - j2; i2++ ) { for ( j = 0; j < lp2h; j++ ) { for ( i = 0; i < lp2h; i++ ) { tmtoe[i+j*lp2h] = tmtoe[i+j*lp2h] + to2[i2+i*lp2h] * mom[j2+i2*lp2h] * te1[j2+j*lp2h]; } } } } tmteo = ( double * ) malloc ( lp3h * lp1h * sizeof ( double ) ); for ( j = 0; j < lp1h; j++ ) { for ( i = 0; i < lp3h; i++ ) { tmteo[i+j*lp3h] = 0.0; } } for ( j2 = 0; j2 < lp2h; j2++ ) { for ( i2 = 0; i2 < lp2h - j2; i2++ ) { for ( j = 0; j < lp1h; j++ ) { for ( i = 0; i < lp3h; i++ ) { tmteo[i+j*lp3h] = tmteo[i+j*lp3h] + te2[i2+i*lp2h] * mom[j2+i2*lp2h] * to1[j2+j*lp2h]; } } } } /* W1 and W2: Interpolation weight matrices. */ w1 = ( double * ) malloc ( lp2h * lp2h * sizeof ( double ) ); for ( j = 0; j < lp2h; j++ ) { for ( i = 0; i < lp2h; i++ ) { w1[i+j*lp2h] = 2.0 / ( double ) ( l * ( l + 1 ) ); } } j = 0; for ( i = 0; i < lp2h; i++ ) { w1[i+j*lp2h] = w1[i+j*lp2h] / 2.0; } if ( ( l % 2 ) == 0 ) { j = lp2h - 1; for ( i = 0; i < lp2h; i++ ) { w1[i+j*lp2h] = w1[i+j*lp2h] / 2.0; } i = lp2h-1; for ( j = 0; j < lp2h; j++ ) { w1[i+j*lp2h] = w1[i+j*lp2h] / 2.0; } } w2 = ( double * ) malloc ( lp3h * lp1h * sizeof ( double ) ); for ( j = 0; j < lp1h; j++ ) { for ( i = 0; i < lp3h; i++ ) { w2[i+j*lp3h] = 2.0 / ( double ) ( l * ( l + 1 ) ); } } i = 0; for ( j = 0; j < lp1h; j++ ) { w2[i+j*lp3h] = w2[i+j*lp3h] / 2.0; } if ( ( l % 2 ) == 1 ) { i = lp3h - 1; for ( j = 0; j < lp1h; j++ ) { w2[i+j*lp3h] = w2[i+j*lp3h] / 2.0; } j = lp1h - 1; for ( i = 0; i < lp3h; i++ ) { w2[i+j*lp3h] = w2[i+j*lp3h] / 2.0; } } /* Cubature weights as matrices on the subgrids. */ for ( j = 0; j < lp2h; j++ ) { for ( i = 0; i < lp2h; i++ ) { w1[i+j*lp2h] = w1[i+j*lp2h] * tmtoe[i+j*lp2h]; } } for ( j = 0; j < lp1h; j++ ) { for ( i = 0; i < lp3h; i++ ) { w2[i+j*lp3h] = w2[i+j*lp3h] * tmteo[i+j*lp3h]; } } /* Pack weight matrices W1 and W2 into the vector W. */ if ( ( l % 2 ) == 0 ) { for ( j = 0; j < lp2h; j++ ) { for ( i = 0; i < lp2h; i++ ) { w[i+2*j*lp2h] = w1[i+j*lp2h]; } } for ( j = 0; j < lp1h; j++ ) { for ( i = 0; i < lp3h; i++ ) { w[i+(2*j+1)*lp2h] = w2[i+j*lp3h]; } } } else { for ( j = 0; j < lp1h; j++ ) { for ( i = 0; i < lp2h; i++ ) { w[i+j*(l+2)] = w1[i+j*lp2h]; } } for ( j = 0; j < lp1h; j++ ) { for ( i = 0; i < lp3h; i++ ) { w[i+lp2h+j*(l+2)] = w2[i+j*lp3h]; } } } /* Free memory. */ free ( te1 ); free ( te2 ); free ( tmteo ); free ( tmtoe ); free ( to1 ); free ( to2 ); free ( w1 ); free ( w2 ); return w; } /******************************************************************************/ void r8mat_transpose_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, char *TITLE, a title. */ { r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 10 September 2013 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, int ILO, JLO, the first row and column to print. Input, int IHI, JHI, the last row and column to print. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2; int i2hi; int i2lo; int i2lo_hi; int i2lo_lo; int inc; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } if ( ilo < 1 ) { i2lo_lo = 1; } else { i2lo_lo = ilo; } if ( ihi < m ) { i2lo_hi = m; } else { i2lo_hi = ihi; } for ( i2lo = i2lo_lo; i2lo <= i2lo_hi; i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; if ( m < i2hi ) { i2hi = m; } if ( ihi < i2hi ) { i2hi = ihi; } inc = i2hi + 1 - i2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, " Row:" ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, " %7d ", i - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Col\n" ); fprintf ( stdout, "\n" ); if ( jlo < 1 ) { j2lo = 1; } else { j2lo = jlo; } if ( n < jhi ) { j2hi = n; } else { j2hi = jhi; } for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, "%5d:", j - 1 ); for ( i2 = 1; i2 <= inc; i2++ ) { i = i2lo - 1 + i2; fprintf ( stdout, " %14g", a[(i-1)+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ double *r8vec_copy_new ( int n, double a1[] ) /******************************************************************************/ /* Purpose: R8VEC_COPY_NEW copies an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, double A1[N], the vector to be copied. Output, double R8VEC_COPY_NEW[N], the copy of A1. */ { double *a2; int i; a2 = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return a2; } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_PRINT prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14g\n", i, a[i] ); } return; } /******************************************************************************/ void r8vec_reverse ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_REVERSE reverses the elements of an R8VEC. Discussion: An R8VEC is a vector of R8's. Example: Input: N = 5, A = ( 11.0, 12.0, 13.0, 14.0, 15.0 ). Output: A = ( 15.0, 14.0, 13.0, 12.0, 11.0 ). Licensing: This code is distributed under the MIT license. Modified: 23 August 2010 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input/output, double A[N], the array to be reversed. */ { int i; int i_hi; double temp; i_hi = n / 2; for ( i = 1; i <= i_hi; i++ ) { temp = a[i-1]; a[i-1] = a[n-i]; a[n-i] = temp; } return; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }