/* Routines in C for general or very specific use: */





In the link above you will find this steps:


Fortran and C interoperability

(a.k.a. How to call C from Fortran, or to call Fortran from C?)

To call C functions from a Fortran program or to call Fortran subroutines from a C program usually involves resolving the following issues:
  • case sensitivity (since C is case sensitive but Fortran is not)
  • possible additional trailing underscore ("_") in C function names (the underscore is added by some Fortran compilers)
  • possible additional Fortran libraries that need to be linked to build the executable
  • Fortran subroutines are always "call by reference" while C can be either "call by reference" or "call by value".In the followings, specific examples are given for Intel compilers and GNU compilers. If you are using other compilers, some changes may be required.
    If you still encounter problems after following the examples below, please email atshpc@ucla.edu

    Fortran calling C functions

Example programs:

Fortran main program, in free format (f-call-c.f):
> program f_call_c
> integer :: i=1, ierr
> double precision :: x=3.14159
> print *, "Fortran calling C, passing"
> print *, "i=",i,"x=",x
> ierr = cfun(i,x)
> end program f_call_c
> 
> [[code]]
C function ([[http://www.ats.ucla.edu/clusters/common/computing/c-fortran/cfun.c|cfun.c]]) -- note the trailing underscore in the function name cfun_:
#include <stdio.h>
int cfun_(int *ip, double *xp)
{
int i = *ip;
double x = *xp;
printf("This is in C function...\n");
printf("i = %d, x = %g\n", i, x);
return 0;
}
> ===How to compile=== 
====Using Intel compilers: icc and ifort:==== 
  • $ ifort -free -c f-call-c.f
  • $ icc -c cfun.c
  • $ ifort -o f-call-c-icc f-call-c.o cfun.o
  • $ ./f-call-c-icc

  • ====Using GNU compilers: gcc and gfortran==== 
$ gfortran -ffree-form -c f-call-c.f
$ gcc -c cfun.c
$ gfortran -o f-call-c-gcc f-call-c.o cfun.o
$ ./f-call-c-gcc
> ==C calling Fortran suboutines== 
===Example programs=== 
C main program ([[http://www.ats.ucla.edu/clusters/common/computing/c-fortran/c-call-f.c|c-call-f.c]]) -- note the trailing underscore in the function name for_sub_:
  • #include <stdio.h>
  • void for_sub_(int *, double *, double *, int *);
  • int main()
  • {
  • int i=1, ny=3;
  • double x=3.14159;
  • double y[]={1.1, 2.2, 3.3};
  • printf("C calling Fortran subroutine, passing\n");
  • printf("i=%d, x=%g\n", i,x);
  • printf("y[]=%g,%g,%g\n",y[0],y[1],y[2]);
  • for_sub_(&i, &x, y, &ny);
  • return 0;
  • }

  • Fortran suboutine, in free format ( [[http://www.ats.ucla.edu/clusters/common/computing/c-fortran/f_sub.f|f_sub.f]]):
subroutine FOR_SUB(i,x,y,ny)
implicit none
integer :: i, ny
double precision :: x
double precision :: y(ny)
print *, "This is in Fortran routine..."
print *, "i = ", i, ", x = ", x
print *, "y = ", y(1:ny)
end subroutine FOR_SUB
> ===How to compile=== 
====Using Intel compilers: icc and ifort==== 
  • $ icc -c c-call-f.c
  • $ ifort -free -c f_sub.f
  • $ icc -o c-call-f-icc c-call-f.o f_sub.o -L/u/local/compilers/intel/fce/current/lib -lifcore -limf
  • $ ./c-call-f-icc

  • ====Using GNU compilers: gcc and gfortran==== 
$ gcc -c c-call-f.c
$ gfortran -ffree-form -c f_sub.f
$ gcc -o c-call-f-gcc c-call-f.o f_sub.o -lgfortran
$ ./c-call-f-gcc
</span>
 
----
** Define the analogous of type "logical" true/false used in Fortran into C (C does not have it defined in the standard library):
 
typedef enum tag_bool
{
FALSE = 0,
TRUE = 1
} logical;
 
----
 
** **CountNRows.c**
** [[http://mywiki-science.wikispaces.com/file/view/CountNRows.c|CountNRows.c]]
** This function reads in a numeric file of the type shown below, and reads ONLY the number of rows
** that contain a numeric value. If the file ends with a set of blank rows, these last rows are not counted.
1 5 2.0
5 6 75.0
8.2 5 3
blank line
blank line
...
 
[[code format="c"]]
//*************************************************************************************************
    int CountNRows(char *input_file)
// This function counts the number of rows of a file and stops at the last numeric value.
// If your file terminates with more than one blank line, this function does NOT
// count those additional lines.
// simone.marras@gmail.com
//
//*************************************************************************************************
{
    int nrows;
    int  i, j, cntr, line_cntr;
    char line[256];
    FILE *file_id;
    /********************************************************************************/
    //Open file and read it throughout once:
    /*FIRST opening of the file to read the number of elements for each category*/
    /* File *. */
    if((file_ID = fopen(input_file, "r")) == NULL){
        printf(" ERROR in CountNRows ! The file %s could not be open\n");
        printf(" The program will EXIT now.\n\n", input_file);
        exit(1);
    }
    else{
        rewind(file_ID); //This line resets the pointer to the beginning of the file
        line_cntr = 0;
        while(fgets(line, sizeof(line), file_ID)!= NULL) {
            if(line[0] < '0')  //zero, not Oh: This control helps to avoid to count any eventual empty line at the end.
                break;
            line_cntr++;
            // ELEMENTS
                nrows = line_cntr ;
            //End else
        }//End WHILE
    }
    fclose(file_ID);
return nrows;
} //END CountNRows()


    • UMFPACK: Full Example of use of the function umfpack_di_triple_to_col.
    • Author of the code: this is taken from the page of Dr. Rouben M. A. Rostamian from the Uni. of Maryland Baltimore.

/* Sample program for the UMFPACK sparse matrix solver
   This is just like the example in UMFPACK's User Guide but we create
   the matrix in the "triplet form" then call umfpack_di_triplet_to_col()
   to convert to "compressed-column form".
   The triplet form is suitable for assembling the mass matrix in
   finite elements.
   Solves the system Ax=b, where:
   A =
       2   3   0   0   0
       3   0   4   0   6
       0  -1  -3   2   0
       0   0   1   0   0
       0   4   2   0   1
   b = (8, 45, -3, 3, 19)
   The solution x is:
   x = (1, 2, 3, 4, 5)
   RR, December 2003
*/
#include <stdio.h>
#include <stdlib.h>
#include <umfpack/umfpack.h>
#include "matvec.h"
#define OK        0
#define FAILED        1
/* Creates a matrix in the triplet form.  Converts it into the
   compressed-column form and stores it in Ap, Ai, Ax which must
   be allocated before calling this function.
   On failure prints a diagnostic to stderr and returns FAILED.
   Otherwise returns OK.
*/
static int create_system(int *Ap, int *Ai, double *Ax, double *b,
        size_t n, size_t nz)
{
    int *Ti, *Tj;
    double *Tx;
    size_t i;
    int status;
    Ti = ivector(nz);
    Tj = ivector(nz);
    Tx = dvector(nz);
    if (Ti == NULL || Tj == NULL || Tx == NULL) {
        fprintf(stderr, "out of memory in function "
                "%s(), file %s, line %d\n",
                __func__, __FILE__, __LINE__);
        status = FAILED;
        goto cleanup;
    }
    i = 0;
    Ti[i] = 0; Tj[i] = 0; Tx[i] = 2.0;
    i++;
    Ti[i] = 1; Tj[i] = 0; Tx[i] = 3.0;
    i++;
    Ti[i] = 0; Tj[i] = 1; Tx[i] = 3.0;
    i++;
    Ti[i] = 2; Tj[i] = 1; Tx[i] = -1.0;
    i++;
    /* we split the (4,1) entry into 1.0+3.0 for illustration */
    Ti[i] = 4; Tj[i] = 1; Tx[i] = 1.0;
    i++;
    Ti[i] = 1; Tj[i] = 2; Tx[i] = 4.0;
    i++;
    Ti[i] = 2; Tj[i] = 2; Tx[i] = -3.0;
    i++;
    Ti[i] = 3; Tj[i] = 2; Tx[i] = 1.0;
    i++;
    Ti[i] = 4; Tj[i] = 2; Tx[i] = 2.0;
    i++;
    Ti[i] = 2; Tj[i] = 3; Tx[i] = 2.0;
    i++;
    Ti[i] = 1; Tj[i] = 4; Tx[i] = 6.0;
    i++;
    Ti[i] = 4; Tj[i] = 4; Tx[i] = 1.0;
    i++;
    /* the second entry for (4,1) */
    Ti[i] = 4; Tj[i] = 1; Tx[i] = 3.0;
    i++;
    /* the right-hand side vector in Ax=b */
    i = 0;
    b[i++] = 8.0;
    b[i++] = 45.0;
    b[i++] = -3.0;
    b[i++] = 3.0;
    b[i++] = 19.0;
    /* convert matrix from triplet form to compressed-column form */
    status = umfpack_di_triplet_to_col(n, n, nz, Ti, Tj, Tx,
            Ap, Ai, Ax, NULL);
    if (status != UMFPACK_OK) {
        fprintf(stderr, "in file %s, line %d: the call to "
                "umfpack_di_triplet_to_col() failed\n",
                __FILE__, __LINE__);
        status = FAILED;
        goto cleanup;
    }
    status = OK;
cleanup:
    free_ivector(Ti);
    free_ivector(Tj);
    free_dvector(Tx);
    return status;
}
int main(void)
{
    void *Symbolic = NULL;
    void *Numeric = NULL;
    int *Ap, *Ai;
    double *Ax;
    double *x, *b;    /* Ax=b */
    size_t n = 5;            /* matrix is nxn */
    size_t nz = 13;            /* length of vectors Ti, Tj, Tx */
    size_t i;
    int status;
    Ap = ivector(n+1);
    Ai = ivector(nz);
    Ax = dvector(nz);
    /* vectors for Ax=b */
    x = dvector(n);
    b = dvector(n);
    if (Ap == NULL || Ai == NULL || Ax == NULL
            || x == NULL || b == NULL) {
        fprintf(stderr, "out of memory in function "
                "%s(), file %s, line %d\n",
                __func__, __FILE__, __LINE__);
        status = EXIT_FAILURE;
        goto cleanup;
    }
    status = create_system(Ap, Ai, Ax, b, n, nz);
    if (status != OK) {
        status = EXIT_FAILURE;
        goto cleanup;    /* diagnostics already printed */
    }
    /* symbolic analysis */
    status = umfpack_di_symbolic(n, n, Ap, Ai, Ax, &Symbolic, NULL, NULL);
    if (status != UMFPACK_OK) {
        fprintf(stderr, "in file %s, line %d: the call to "
                "umfpack_di_symbolic() failed\n",
                __FILE__, __LINE__);
        status = EXIT_FAILURE;
        goto cleanup;
    }
    /* LU factorization */
    status = umfpack_di_numeric(Ap, Ai, Ax, Symbolic,
            &Numeric, NULL, NULL);
    if (status != UMFPACK_OK) {
        if (status == UMFPACK_WARNING_singular_matrix)
            fprintf(stderr, "according to umfpack_di_numeric() "
                    "matrix is singluar.\n"
                    "computation aborted\n");
        else
            fprintf(stderr, "in file %s, line %d: the call to "
                    "umfpack_di_numeric() failed\n",
                    __FILE__, __LINE__);
        status = EXIT_FAILURE;
        goto cleanup;
    }
    /* solve system */
    status = umfpack_di_solve(UMFPACK_A, Ap, Ai, Ax, x, b,
            Numeric, NULL, NULL);
    if (status != UMFPACK_OK) {
        fprintf(stderr, "in file %s, line %d: the call to "
                "umfpack_di_solve() failed\n",
                __FILE__, __LINE__);
        status = EXIT_FAILURE;
        goto cleanup;
    }
    for (i = 0; i < n; i++)
        printf("x[%d] = %g\n", i, x[i]);
    status = EXIT_SUCCESS;
cleanup:
    umfpack_di_free_numeric(&Numeric);
    umfpack_di_free_symbolic(&Symbolic);
    free_ivector(Ap);
    free_ivector(Ai);
    free_dvector(Ax);
    free_dvector(x);
    free_dvector(b);
    return status;
}



    • Small library for dynamic allocation with the option of standard C 0-indexing, or of1-indexing (as in Fortran)

    1. Source code: DYNALLOC.c DYNALLOC.h

Example:
//Integer matrix
int **iarray2d(size_t nrows, size_t ncols, int index)
{
    //nrows = number of rows
    //ncols = ncolumns
    //index = 1 if you want a fortran-type indexing from 1
    //index = 0 if you want the standard C allocation starting from 0
    long i;
    int **mat;
    //Malloc pointers to rows
    mat= malloc((nrows+index+1)*sizeof(int*));
    if (!mat){
        printf("Matrix MAT could not be allocted");
        exit(1);
    }
    mat[0] = NULL;
    for(i=index; i<nrows+index+1; i++){
        mat[i] = malloc((ncols+index)*sizeof(int*));
        if (!mat[i]){
            printf("Matrix MAT could not be allocted\n The exec. will Exit now\n");
            exit(1);
        }
    }
return mat;
}
//FREE MEMORY:
//Integer matrix
void free_iarray2d(int **mat, size_t nrows, size_t ncols, int index)
{
    //nrows = number of rows
    //ncols = ncolumns
    //index = 1 if you want allocation in the fortran way starting from 1
    //index = 0 if you want the standard C allocation way starting from 0
    long i;
    for(i=index; i<nrows+index+1; i++)
        free(mat[i]);
    free(mat);
return;
}


    • STRUCTS: Dynamic allocation of the arguments (arrays) of a generic struct

    1. READ_BC.h
typedef struct boundary_conditions
{
    unsigned int num_nodes; // Number of boundary nodes
    unsigned int *node;     // Node index number
    unsigned int *code;     // Array of the boundary codes
    float *values;          // Array of boundary values: size [nbdy_nodes]x[nvariables]
} boundary_conditions;
void read_boundary_conditions(char *bdy_file, boundary_conditions *bdy_conds);
    1. READ_BC.c (Dynamic allocation is done here)
void READ_BC(char *bdy_file, boundary_conditions *BDYCS)
{
    FILE *bdyf_ID;
    char header[16];
    bdy_conds->node = malloc(10 * sizeof(*bdy_conds->node));
    free(BDYCS.NODE);
return;
}


    • Function that counts the number of points surrounding a point of a computational mesh

This function comes from the algorithm of Lohner explained in his book [1]. I report its C implementation here (for a matlab version go here ), with a couple of explanations on its use and the type of mesh that it accepts.

[1] Lohner, Rainald (2008),'Applied computational fluid dynamics techniques', John Wiley and Sons
Source: PSUP.c and PSUP.h
You also NEED nrutil.c and nrutil.h in your compilation directory
/******************************************************************
*    *PSUP(int **CONN, int conn_flag, int nnodes, int nelem, int max_nnodes, unsigned int array_numb)
*
*    Points Surrounding Point (PSUP):
*
*    Input:
*    - **CONN is a 2D array containing the mesh. The mesh should have been
*      already read within the main file that calls PSUP, and should be stored into CONN.
*      Now, CONN can be of two different forms; see below:
*
*    - CONN: array of the connectivity matrix:
*                        WARNING: this function works ONLY for unstructured 2D grids.
*                        It must contain the conn matrix in the format:
*                        conn(1:nelem, 1:elnnode+conn_flag)
*                        where nelem is the mesh number of elements and
*                        elnnode is the max number of nnodes that any element may have.
*                        conn_flag = 1 means that the 1st column of the CONN matrix represents
*                                      the element number (see below);
*                        conn_flag = 0 means that the 1st column of the CONN matrix is already
*                                      a node number (see below);
*                        ex with conn_flag = 1:
*                        1 2 5 3
*                        2 6 3 1
*                        3 3 1 8
*                        4 9 8 5
*                        5 ... ... ...
*                        ...
*                        ex with conn_flag = 0:
*                        2 5 3
*                        6 3 1
*                        3 1 8
*                        9 8 5
*                        ...
*
*
*    - nnodes:     Number of points in the mesh
*    - nelem:      Number of elements in the mesh
*    - max_nnodes: Maximum number of nodes of the elements: ex. if a mesh is made of quadrilaterals,
*                         max_nnodes = 4; if made of triangles, max_nnodes = 3.
*    - array_numb: use 1 of you allocated CONN on a 1-base index (as in Fortran),
*                  or 0 if you allocated on a 0-index base (standard C)
*
*    Output (returns):
*    - An array of boundary flags: this array contains the nudes that are on the mountain boundary.
*        WARNING: This array is correct ONLY if the mesh in input is made of quadrilaterals.
*        For this to work on triangles as well, an additional algorithm must be added to this
*        function; the additional algorithm is also found in Lohner [1].
*
*    Revised and fixed: June 20 2009. Working correctly.
*    Simone Marras
*
*    REFERENCES:
*    [1] Lohner,R. (2008), Applied computational fluid dynamics techniques",ed.II, John Wiley and Sons
*
******************************************************************/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include "nrutil.h"
int *PSUP(int **CONN, int conn_flag, int nnodes, int nelem, int max_nnodes, unsigned int array_numb)
{
    int bnodes, *BFLAG;
    int i,j, ipoin, jpoin, inode, ielem, ipoi1, ipoi2, istor;
    int szesup1;
    int iesup;
    int max_psup;
    int *esup1, *esup2, *nesup, *psup1, *psup2, *npsup;
    int *lpoin;
    int **CONN_local;
    int TEMP;
    //Dynamic allocation
    CONN_local = imatrix(0,nelem+1, 0,max_nnodes+1);
    esup2 = calloc(nnodes+1,sizeof(int*));
    esup1 = calloc(1,sizeof(int*));    // Will be reallocated only later
    nesup = calloc(nnodes+1,sizeof(int*));
    psup2 = calloc(nnodes+1,sizeof(int*));
    lpoin = calloc(nnodes+1,sizeof(int*));
    psup1 = calloc(1,sizeof(int*)); // Will be reallocated only later
    npsup = calloc(nnodes+1,sizeof(int*));
    BFLAG = calloc(nnodes+1,sizeof(int *)); // Will be reallocated only later
    //IMPORTANT NOTICE: BFLAG is 0-indexed BUT the valuable indexes for node flagging start with 1!!!
    //                    BFLAG[0] has NO MEANING at all.
    if(esup1==NULL || esup2==NULL || nesup==NULL || psup1==NULL || psup2==NULL || npsup==NULL || lpoin==NULL) {
    printf(" # !!! ERROR IN DYNAMIC ALLOCATION\n THE PROGRAM WILL EXIT NOW! ");
    exit(1);
    }
    printf("\n###############################\n# Elements surrounding points #\n###############################\n");
    printf("\n # nodes %d, nelem %d, max_nnodes %d\n", nnodes, nelem, max_nnodes);
    //Pre-processing of the connectivity matrix. Done to set the matrix into the right shape used next:
    if( conn_flag == 1){
        for(i=array_numb; i<=nelem-1+array_numb; i++){
            for(j=array_numb; j<=max_nnodes+array_numb; j++){
                TEMP = CONN[i][j+1];
                CONN_local[i][j] = TEMP;
                //printf("%d ", CONN_local[i][j]);
            }
            //printf("\n");
        }
    }
    else{
        for(i=array_numb; i<=nelem-1+array_numb; i++){
            for(j=array_numb; j<=max_nnodes+array_numb; j++){
                CONN_local[i][j] = CONN[i][j];
            }
        }
    }
/* LOHNER's ALGORITHM BEGINS HERE */
//Element pass 1: count the number of elements connected to each point:
    for(ielem=array_numb; ielem<nelem+array_numb; ielem++){
        for(inode=array_numb; inode<max_nnodes+array_numb; inode++){
        //Update storage counter, storing ahead:
            ipoi1 = CONN_local[ielem][inode] + 1;
            esup2[ipoi1-1+array_numb] = esup2[ipoi1-1+array_numb] + 1;
            //printf("PASS 1 esup2[%d] = %d\n", ipoi1, esup2[ipoi1]); //OK until here
        }
    }
//Storage/reshugffling pass 1:
    for(ipoin=array_numb+1; ipoin<nnodes+1+array_numb; ipoin++){
    //Update storage counter and store:
        esup2[ipoin-1+array_numb] = esup2[ipoin-1+array_numb] + esup2[ipoin-1+array_numb-1]; //Ok
        //printf("esup2[%d] = %d\n", ipoin, esup2[ipoin]); //Ok
    }
    szesup1 = esup2[nnodes+array_numb];
    esup1 = realloc(esup1,(szesup1+1)*sizeof(int*));
    //Set esup1 to zero:
    for(i=array_numb; i<szesup1+1; i++)
        esup1[i] = 0;
    if(esup1 != NULL)
        printf(" #\n # esup1 Now reallocating more memory... \n");
    else
        printf(" !!! NOT ENOUGH MEMORY \n");
//Element pass 2: store the elements in esup1:
    for(ielem=array_numb; ielem<nelem+array_numb; ielem++){
        for(inode=array_numb; inode<max_nnodes+array_numb; inode++){
        //Update storage counter, storing in esup1:
            ipoin = CONN_local[ielem][inode]; //Ok
//            printf("ipoin = %d\n", ipoin);
            istor = esup2[ipoin-1+array_numb] + 1; //Ok
//            printf("istor = %d\n", istor);
            esup2[ipoin-1+array_numb] = istor+1-array_numb; //OK
            esup1[istor-1+array_numb] = ielem+1-array_numb; //Ok
//            printf("esup2[%d] = %d\n", ipoin-1+array_numb, esup2[ipoin-1+array_numb]); //Ok
//            printf("esup1[%d] = %d\n", istor-1+array_numb, esup1[istor-1+array_numb]); //Ok
        }
    }
//Storage/reshuffling pass 2 and count the number of elements sorrounding point ipoin:
    for(ipoin = nnodes+array_numb; ipoin>=array_numb+1; ipoin--){
        nesup[ipoin] = esup2[ipoin] - esup2[ipoin-1];
        esup2[ipoin] = esup2[ipoin-1]; //OK
    }
        esup2[array_numb] = 0; //OK
        nesup[array_numb] = esup2[array_numb+1] - esup2[array_numb]; //OK
    //for(i=array_numb; i<=nnodes+array_numb; i++){
    //    printf("nesup2[%d] = %d\n", i, nesup[i]);//OK
    //    printf("esup2[%d] = %d\n", i, esup2[i]);//OK
    //}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Points surrounding points:
printf("\n ###############################\n # Points surrounding points #\n");
    psup2[array_numb] = 0;
    istor = 0;
    //!!! IMPORTANT NOTE on the ALGORITHM: Here you will find a set of TWO loops practically identical:
    //the first time that the loop is used, it obtains the number of elements that
    //will fill the array "psup1" that until now was not yet declared and allocated;
    //The second time, it actually fills in such array (this does not happen in the
    //original algorithm of Lohner presented in the book.):
    //1) FIRST OF 2 LONG LOOPS to obtain the size of "psup1"
    for(ipoin=array_numb; ipoin<nnodes+array_numb; ipoin++){
    //Loop over the elements sorrounding the point:
        for(iesup=esup2[ipoin]+1; iesup<=esup2[ipoin+1]; iesup++){
            ielem = esup1[iesup-1+array_numb]; //OK
            //printf("iesup = %d  ielem = %d\n", iesup, ielem);//OK
            //Loop over the nodes of the element
            for(inode=array_numb; inode<max_nnodes+array_numb; inode++){
                jpoin = CONN_local[ielem-1+array_numb][inode-1+array_numb];
                //printf("jpoin[%d] = %d\n", inode, jpoin); //OK
                if(jpoin != ipoin && lpoin[jpoin-1+array_numb] != ipoin){//OK
                    //Update storage counter, storing ahead, and mark lpoin:
                    istor = istor + 1;
                    //printf(" ISTOR = %d\n ", istor);//OK
                    //NOTE: The original loop of the algorithm (corresponding to our 2nd loop below this),
                    //here already fills in psup1! We do it in the next after allocating it.
                    lpoin[jpoin-1+array_numb] = ipoin; //OK
                    //printf(" lpoin[%d] = %d\n ",jpoin-1+array_numb, lpoin[jpoin-1+array_numb]);
                }
            }
        }
    }
    //Reallocate psup1 with istor elements and BFLAG with n bnodes:
    psup1 = realloc(psup1,(istor+1)*sizeof(int*));        //OK
    lpoin = realloc(lpoin, (nnodes+1)*sizeof(int*));    //OK
    if(psup1!=NULL)
        printf(" #\n # psup1 Now reallocating more memory... \n");
    else
        printf("\n !!! NOT ENOUGH MEMORY to allocate psup1\n");
    //Null the used arrays before reusing it:
    //OK
    for(i=0; i<nnodes+1; i++){
        lpoin[i] = 0;
        psup2[i] = 0;
        npsup[i] = 0;
    }
    for(i=0; i<istor+1; i++){
        psup1[i] = 0;
    }
    istor = 0;
    //2) SECOND OF 2 LONG LOOPS to finally FILL IN "psup1"
    for(ipoin=array_numb; ipoin<nnodes+array_numb; ipoin++){
    //Loop over the elements sorrounding the point:
        for(iesup=esup2[ipoin]+1; iesup<=esup2[ipoin+1]; iesup++){
            ielem = esup1[iesup-1+array_numb]; //OK
            //printf("iesup = %d  ielem = %d\n", iesup, ielem);//OK
        //Loop over the nodes of the element
            for(inode=array_numb; inode<max_nnodes+array_numb; inode++){
                jpoin = CONN_local[ielem-1+array_numb][inode-1+array_numb];
                //printf("jpoin[%d] = %d\n", inode, jpoin); //OK
                if(jpoin != ipoin && lpoin[jpoin-1+array_numb] != ipoin){//OK
                    //Update storage counter, storing ahead, and mark lpoin:
                    istor = istor + 1;
                    //printf(" ISTOR = %d\n ", istor);//OK
                    psup1[istor-1+array_numb] = jpoin-1+array_numb;//OK
                    //printf(" psup1[%d] = %d\n ",istor-1+array_numb, psup1[istor-1+array_numb]);//OK
                    lpoin[jpoin-1+array_numb] = ipoin-1+array_numb;//OK
                    //printf(" lpoin[%d] = %d\n ",jpoin-1+array_numb, lpoin[jpoin-1+array_numb]);//OK
                }
            }
        }
        //Update storage counters:
            psup2[ipoin+1] = istor;//OK
            npsup[ipoin] = psup2[ipoin+1]-psup2[ipoin];//OK
            //printf("psup2[%d] = %d\n", ipoin+1,psup2[ipoin+1]); //OK
            printf(" # Node %d is sorrounded by %d nodes\n", ipoin,npsup[ipoin]); //OK
            //Here the Boundary flag is filled with only the boundary nodes.
            //Be very careful because this is correct only for quadrilaterals!!!
            //If you want the b.flag to be filled correctly when triangles are used, you
            //need to add another part of algorithm following Lohner: chp. 2.2.3, pag.12
            // 'elements sorrounding elements'
    }//END second loop
    //Get the maximum number of points surrounding points:
    // This can be used as the maximum number of nonzero elements
    // appearing among all the rows of the stiffness matrix
        max_psup = 0;
        for(ipoin=0; ipoin<nnodes+1; ipoin++){
            if(npsup[ipoin] > max_psup){
                //printf(" npsup[%d] = %d\n", ipoin, npsup[ipoin]);//OK
                max_psup = npsup[ipoin];//OK
            }
        }
        printf(" # max_psup = %d\n ", max_psup);
        //Here we count the number of nodes that are boundary nodes. This count is
        //necessary to allocate the memory for the Boundary FLAG.
        //IMPORTANT NOTICE: BFLAG is 0-indexed BUT the valuable indexes for node flagging start with 1!!!
        //                    BFLAG[0] has NO MEANING at all.
        for(ipoin = 0; ipoin<nnodes+1; ipoin++){
            if( max_nnodes == 3 && npsup[ipoin] < 6)
            //Boundary node for TRI element grids
            {
                BFLAG[ipoin] = 0;
            }
            else if( max_nnodes == 4 &&  npsup[ipoin] < 8)
            //Boundary node for QUAD element grids
            {
                BFLAG[ipoin] = 0;
            }
            else //if internal node
                BFLAG[ipoin] = npsup[ipoin];
            printf(" BFLAG[%d] = %d\n", ipoin, BFLAG[ipoin]);
        }
printf(" #\n # Maximum number of points surrounding points: %d\n###############################\n\n", max_psup);
/*Free memory*/
    free_imatrix(CONN_local, 0,nelem+1, 0,max_nnodes+1);
    free(esup2);
    free(esup1);
    free(nesup);
    free(psup2);
    free(lpoin);
    free(psup1);
    free(npsup);
    //free(BFLAG);    // NOTICE: Uncomment this if you are not returning BFLAG to the calling function. Otherwise, keep it commented
                    // but DO NOT FORGET to free BFLAG in the calling function after use of BFLAG.
return BFLAG;
}



    1. Sorting in ascending order 1D and 2D arrays (of floats): (go here for its Matlab and Octave version)
    • 1D array of size 'rows':
    • void fsortAscend(float *array, int length, int array_numb)
      {
        int i, j, temp, test;
       
          for(i = length + array_numb - 1; i > 0+array_numb; i--)
          {
          test=0;
              for(j = 0; j < i; j++)
              {
                if(array[j] > array[j+1]) /* compare neighboring elements */
                {
                  temp = array[j];    /* swap array[j] and array[j+1] */
                  array[j] = array[j+1];
                  array[j+1] = temp;
                  test=1;
                }
              } /*end for j*/
              if(test==0) break; /*will exit if the list is sorted!*/
          } /*end for i*/
       }
    • 2D array of size 'rows' by 'cols': this function sorts a 2D array row by row; column elements are not sorted in any way; see the ex. below:
    • Before:
      300  4 6 7
      200 1 4 2
      10   5 7 32
      150 6 2 4
      After:
      10   5 7 32
      150 6 2 4
      200 1 4 2
      300  4 6 7
    • void fsortAscend2D(float **array, int rows, int cols, int array_numb)
      {
        int i,j,k, test;
       
        float *temp2d;
        float *temp;
       
        temp = malloc((cols+1)*sizeof(float*));
       
          for(i = rows + array_numb - 1; i > 0+array_numb; i--)
          {
          test=0;
              for(j = 1 ; j < i; j++)
              {
                  if(array[j][1] > array[j+1][1]) /* compare neighboring elements */
                  {
                      for(k=array_numb; k<cols+array_numb; k++){
                          temp[k] = array[j][k];    /* swap array[j] and array[j+1] */
                          array[j][k] = array[j+1][k];
                          array[j+1][k] = temp[k];
                          test=1;
                      }
                  }
              } /*end for j*/
              if(test==0) break; /*will exit if the list is sorted!*/
          } /*end for i*/
       
          free(temp);
      }
See also the cboard forum for more: **cboard Forum**


    • Dynamic memory allocation: dynamic allocation of arrays is quite straight forward and you can find it anywhere. Here I want to add how to allocate dynamically an array of pointers to chars; nominally, something that, if allocated statically, would look like this:

const char *mystrings[4] = {"hello", "how", "are", "you?"};
where 4 is the number of strings that *mystring would contiain, and the size of each of them is given by each entries.

However, more than likely you may need to allocate it dynamically both in the number of strings that mystrings[] contains, and in terms of the size of each entry string as well:

Here the way I did it: it seems to be working correctly, but if you find a mistake that may not be visible until debugging, please, login and correct it

/* Declarations */
#define MAX_STR_LENGTH 24
char **strings;
char str_name[MAX_STR_LENGTH];
int i, n_strings;
/*Miscellaneous Declarations and assignments*/
    n_strings = whatever number;
    strcpy(str_name, "hello");
/* Dynamic allocation */
     strings = (char**) malloc(n_strings* sizeof(char *));
    for(i = 0;  i < nstrings; i++){
          strings[i] = (char*) malloc( sizeof(str_name) * sizeof(char *));
    }
/* Now we can assign each entry with a string name as follows */
     for(i = 0;  i < nstrings; i++){
          strncpy(strings[i], str_name, sizeof(str_name));
     }
/* FREE allocated memory */
   free(strings);