Re: Segmentation fault when using the fortran subroutine

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Thu Oct 24 2013 - 21:45:39 MDT

Not sure of the problem. Your wrapped code and the NCL script look good.
Some assorted comments:

[1]

The fortran documentation states:
   F must be dimensioned at least (M+1)*(N+1).

Are the arguments the correct sizes?

---
[2]
I know that NCL is distributed with LAPACK. However,
I did not know NCL was distributed with FISHPACK.
The LAPACK is double precision so if NCL has FISHPACK, it is
likely double precision also. You are passing real values
and not double. This could cause the type of error
you are getting. Please try double precision.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
The following is text extracted from 2 power point slides
discussed in the NCL Workshop
================================================
            Accessing LAPACK  (1 of 2)
  double precision LAPACK (BLAS) => distributed with NCL
       explicitly link LAPACK lib via fortran interface: WRAPIT
  eg: subroutine dgels solves [over/under]determined real linear systems
C NCLFORTSTART
       SUBROUTINE DGELSI( M, N, NRHS, A, B, LWORK, WORK )
       IMPLICIT   NONE
       INTEGER    M, N, NRHS, LWORK
       DOUBLE PRECISION   A( M, N ), B( M, NRHS), WORK(LWORK)
C NCLEND
C                                                     declare local 
variables
       INTEGER     INFO
       CHARACTER*1 TRANS
       TRANS = "N”
       CALL DGELS(TRANS, M,N,NRHS,A,LDA,B,LDB,WORK,LWORK,INFO)
       RETURN
       END
================================================
           Accessing LAPACK  (2 of 2)
external DGELS "./dgels_interface.so”
; NAG example: http://www.nag.com/lapack-ex/node45.html
; These are transposed from the fortran example
  A = (/ (/ -0.57, -1.93, 2.30, -1.93, 0.15, -0.02 /), \  ; (4,6)
           (/ -1.28,  1.08, 0.24,  0.64, 0.30,  1.03 /), \
           (/ -0.39, -0.31, 0.40, -0.66, 0.15, -1.43 /), \
           (/  0.25, -2.14,-0.35,  0.08,-2.13,  0.50 /)  /)*1d0     ; 
must be double
     dimA  = dimsizes(A)
     N     = dimA(0)         ; 4
     M     = dimA(1)        ; 6
     B     = (/-2.67 ,-0.55 ,3.34, -0.77, 0.48, 4.10/)*1d0       ; must 
be double
                                   ; LAPACK wants 2D
     nrhs  = 1
     B2    = conform_dims ( (/nrhs,M/), B, 1 )    ; (1,6)
     B2(0,:) =B
     lwork=500 
; allocate space
     work   = new ( lwork, "double", "No_FillValue")         ; must be 
double
     DGELS::dgelsiw(M, N, nrhs, A , B2, lwork, work )
     print(B2(0,0:N-1))
---
Hope this helps.
On 10/23/13 9:39 AM, Leung YU Ting wrote:
> HI
>
>       When I want to use the HWSSSP of fishpack Fortran subroutine in NCL. I also add a subroutine hstub.f to define the dimension of variables as shown below;
>
> C NCLFORTSTART
>        SUBROUTINE HWSSSPDRV (TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,
>       1                      BDPS,BDPF,ELMBDA,F,NLON,NLAT,PERTRB,
>       2                      IERROR,W,LWORK)
>        DIMENSION       F(NLON,NLAT) ,BDTS(NLON),BDTF(NLON),BDPS(NLAT),
>       1                BDPF(NLAT),W(LWORK)
> C NCLEND
>        print*, "1"
>        CALL HWSSSP(TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS,
>       1            BDPF,ELMBDA,F,NLON,PERTRB,IERROR,W)
>        print*, "2"
>        RETURN
>        print*, "3"
>        END
>
> The ncl script used is shown below;
>
>      rawtime=hdata->time
>      time=ut_calendar(rawtime, -3)
>      lat=hdata->lat(::-1)
>      lon=hdata->lon
>
>      wantedtime=2000120100
>
>      timeind=ind(wantedtime .eq. time)
>
>      wantedhgt1=short2flt(hdata->hgt(timeind-1, {500}, ::-1, :))
>      wantedhgt2=short2flt(hdata->hgt(timeind+1, {500}, ::-1, :))
>
>      guv1=z2geouv(wantedhgt1, lat, lon, 1)
>      gu1=guv1(0, :, :)
>      gv1=guv1(1, :, :)
>      gvort1=uv2vr_cfd(gu1, gv1, lat, lon, 3)
>      copy_VarCoords(wantedhgt1, gvort1)
>
>      guv2=z2geouv(wantedhgt2, lat, lon, 1)
>      gu2=guv2(0, :, :)
>      gv2=guv2(1, :, :)
>      gvort2=uv2vr_cfd(gu2, gv2, lat, lon, 3)
>      copy_VarCoords(wantedhgt2, gvort2)
>
>      gvort_tend=gvort2-gvort1
>      copy_VarCoords(gvort1, gvort_tend)
>
>      PI = 3.14159265358979
>      NLAT = dimsizes(lat)
>      NLON = dimsizes(lon)
>      PERTRB = 0.
>      IERROR = 1
>      TS= 0.
>      TF= PI
>      M=NLAT-1
>      MBDCND=9
>      BDTS=new(NLON, float)
>      BDTF=new(NLON, float)
>      PS=0.
>      PF=2.*PI*357.5/360.
>      N=NLON-1
>      M=NLAT-1
>      NBDCND=0
>      BDPS=new(NLAT, float)
>      BDPF=new(NLAT, float)
>      ELMBDA=0.
>      IDIMF= NLON
>      LWORK = 11*(M+1)+6*(N+1)
>      W=new(LWORK,float)
>      F=gvort_tend
>
>      HSTUB::HWSSSPDRV(TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS,BDPF,ELMBDA,F,NLON,NLAT,PERTRB,IERROR,W,LWORK)
>
> The output is;
>   1
>   2
> Segmentation fault (core dumped)
>
> It seems the script stop at return the variables. So I use gdb to check the script. It is stop at _NclGetObj. What can I do to solve this problem?
>
> Program received signal SIGSEGV, Segmentation fault.
> 0x000000000061d5f2 in _NclGetObj ()
> (gdb) where
> #0  0x000000000061d5f2 in _NclGetObj ()
> #1  0x000000000064828f in VarValRep ()
> #2  0x000000000064de18 in _NclGetVarRepValue ()
> #3  0x0000000000649689 in VarVarWrite ()
> #4  0x000000000064e2d4 in _NclAssignVarToVar ()
> #5  0x00000000006cd4f9 in _NclRemapIntrParameters ()
> #6  0x00000000006c2070 in CallINTRINSIC_PROC_CALL ()
> #7  0x00000000006c91ad in _NclExecute ()
> #8  0x00000000005af1b1 in yyparse ()
> #9  0x00000000005aa592 in main ()
> (gdb)
>
> Regards
> Marco
>   		 	   		
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Oct 24 21:45:49 2013

This archive was generated by hypermail 2.1.8 : Fri Nov 01 2013 - 08:58:14 MDT