i started writing a wrapper for fortran code to calculate svd.
currently, i'm not sure which algorithm would be optimal for my problem.
perhaps svd is actually not the right thing. interpolative decompositions might also work.

in order to be able to load the library into sbcl i had to compile lapack and id_lib with -fPIC as a shared library
this is easily configured for lapack using cmake -i.

in id_lib i added the flags -fPIC -shared to the makefile and in the end executed the following command to create the library:

gfortran -shared -o libid_lib.so  prini.o dfft.o id_rand.o id_rtrans.o idd_snorm.o idd_sfft.o idd_frm.o idd_house.o idd_qrpiv.o idd_id.o idd_svd.o idd_id2svd.o iddr_rid.o iddr_rsvd.o iddp_rid.o iddp_rsvd.o iddr_aid.o iddr_asvd.o iddp_aid.o iddp_asvd.o idz_snorm.o idz_sfft.o idz_frm.o idz_house.o idz_qrpiv.o idz_id.o idz_svd.o idz_id2svd.o idzr_rid.o idzr_rsvd.o idzp_rid.o idzp_rsvd.o idzr_aid.o idzr_asvd.o idzp_aid.o idzp_asvd.o -llapack -L/usr/local/lib

according to the documentation i believe i will need idzr_asvd which depends on the initialization function idzr_aidi

computes the SVD of an arbitrary (generally
dense) matrix, to a specified rank; this routine
is randomized

        integer len
        parameter(len = 1 000 000)

        integer m,n,krank,ier
        real*8 errmax,errrms,s(100 000)
        complex*16 a(len),b(len),u(len),v(len),w(len)

c       Calculate an SVD approximating a.
        call idzr_aidi(m,n,krank,w)
        call idzr_asvd(m,n,a,krank,w,u,v,s,ier)

c       Construct b = u diag(s) v^*.
        call reconsvd(m,krank,u,s,n,v,b)

        subroutine reconsvd(m,krank,u,s,n,v,a)
c
c       forms a = u diag(s) v^*.
c
c       input:
c       m -- first dimension of u and a
c       krank -- size of s, and second dimension of u and v
c       u -- leftmost matrix in the product a = u diag(s) v^*
c       s -- entries on the diagonal in the middle matrix
c            in the product a = u diag(s) v^*
c       n -- second dimension of a and first dimension of v
c       v -- rightmost matrix in the product a = u diag(s) v^*
c
c       output:
c       a -- matrix product u diag(s) v^*
        implicit none
        integer m,n,krank,j,k,l
        real*8 s(krank)
        complex*16 u(m,krank),v(n,krank),a(m,n),sum
        do k = 1,n
          do j = 1,m
            sum = 0
            do l = 1,krank
              sum = sum+u(j,l)*s(l)*conjg(v(k,l))
            enddo ! l
            a(j,k) = sum
          enddo ! j
        enddo ! k
        return
        end

        subroutine idzr_aidi(m,n,krank,w)
c
c       initializes the array w for using routine idzr_aid.
c
c       input:
c       m -- number of rows in the matrix to be ID'd
c       n -- number of columns in the matrix to be ID'd
c       krank -- rank of the ID to be constructed
c
c       output:
c       w -- initialization array for using routine idzr_aid
c
        implicit none
        integer m,n,krank,l,n2
        complex*16 w((2*krank+17)*n+21*m+80)

        subroutine idzr_asvd(m,n,a,krank,w,u,v,s,ier)
c
c       constructs a rank-krank SVD  u diag(s) v^*  approximating a,
c       where u is an m x krank matrix whose columns are orthonormal,
c       v is an n x krank matrix whose columns are orthonormal,
c       and diag(s) is a diagonal krank x krank matrix whose entries
c       are all nonnegative. This routine uses a randomized algorithm.
c
c       input:
c       m -- number of rows in a
c       n -- number of columns in a
c       a -- matrix to be decomposed; the present routine does not
c            alter a
c       krank -- rank of the SVD being constructed
c       w -- initialization array that routine idzr_aidi
c            has constructed (for use in the present routine,
c            w must be at least
c            (2*krank+22)*m+(6*krank+21)*n+8*krank**2+10*krank+90
c            complex*16 elements long)
c
c       output:
c       u -- matrix of orthonormal left singular vectors of a
c       v -- matrix of orthonormal right singular vectors of a
c       s -- array of singular values of a
c       ier -- 0 when the routine terminates successfully;
c              nonzero otherwise
c
c       _N.B._: The algorithm used by this routine is randomized.
c
        implicit none
        integer m,n,krank,lw,ilist,llist,iproj,lproj,icol,lcol,
     1          iwork,lwork,iwinit,lwinit,ier
        real*8 s(krank)
        complex*16 a(m,n),u(m,krank),v(n,krank),
     1             w((2*krank+22)*m+(6*krank+21)*n+8*krank**2
     2              +10*krank+90)
i made a few calls to id_lib's idzr_asvd. it didn't crash but so far i haven't verified that the svd works
6Gb of RAM and 40s solves
m=20000 n=10002 k=300
with random entries in a
185s and 4.5Gb RAM solve
m=10000 n=10000 k=900
3.9Gb RAM and 192s solve
random
m=10000 n=10000 k=1000
the lapack svd seems to run approximately with the same speed as id_dist's svd.
i will have to compare their results
both functions return the same singular values (they differ by 7e-14)