Ignore:
Timestamp:
Jan 28, 2010 6:09:34 PM (13 years ago)
Author:
srkline
Message:

Updated the MonteCarlo? code to allow 4 processors, but simply copying the function 4 times, and defining 4 different random number generators. Still can't figure out what the problem is with threading a single version, but not worth the effort. Copy/paste is way faster.

Also added some simple (non-optimized) calculations for using Debye's sphere method. These are largely undocumented at this point - so see the code. These are XOP versions of the old ipf code I've used in the past, and stripped of the now-obsolete AltiVec? code (I now lose the 4x speedup from the vectorization...)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/XOP_Dev/MonteCarlo/MonteCarlo.c

    r591 r623  
    1111#include "XOPStandardHeaders.h"                 // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h 
    1212#include "MonteCarlo.h" 
     13#include "DebyeSpheres.h" 
    1314 
    1415static int gCallSpinProcess = 1;                // Set to 1 to all user abort (cmd dot) and background processing. 
     
    597598#undef RNMX 
    598599 
     600 
     601//////////a complete copy of ran1(), simply renamed ran1a() 
     602#define IA 16807 
     603#define IM 2147483647 
     604#define AM (1.0/IM) 
     605#define IQ 127773 
     606#define IR 2836 
     607#define NTAB 32 
     608#define NDIV (1+(IM-1)/NTAB) 
     609#define EPS 1.2e-7 
     610#define RNMX (1.0-EPS) 
     611 
     612float ran1a(long *idum) 
     613{ 
     614        int j; 
     615        long k; 
     616        static long iy=0; 
     617        static long iv[NTAB]; 
     618        float temp; 
     619 
     620        if (*idum <= 0 || !iy) { 
     621                if (-(*idum) < 1) *idum=1; 
     622                else *idum = -(*idum); 
     623                for (j=NTAB+7;j>=0;j--) { 
     624                        k=(*idum)/IQ; 
     625                        *idum=IA*(*idum-k*IQ)-IR*k; 
     626                        if (*idum < 0) *idum += IM; 
     627                        if (j < NTAB) iv[j] = *idum; 
     628                } 
     629                iy=iv[0]; 
     630        } 
     631        k=(*idum)/IQ; 
     632        *idum=IA*(*idum-k*IQ)-IR*k; 
     633        if (*idum < 0) *idum += IM; 
     634        j=iy/NDIV; 
     635        iy=iv[j]; 
     636        iv[j] = *idum; 
     637        if ((temp=AM*iy) > RNMX) return RNMX; 
     638        else return temp; 
     639} 
     640#undef IA 
     641#undef IM 
     642#undef AM 
     643#undef IQ 
     644#undef IR 
     645#undef NTAB 
     646#undef NDIV 
     647#undef EPS 
     648#undef RNMX 
     649/////////////// 
     650 
     651 
    599652//////////////////////// 
    600653#define MBIG 1000000000 
     
    645698#undef FAC 
    646699 
     700//////////////////////// a complete copy of ran3() renamed ran3a() 
     701#define MBIG 1000000000 
     702#define MSEED 161803398 
     703#define MZ 0 
     704#define FAC (1.0/MBIG) 
     705 
     706float ran3a(long *idum) 
     707{ 
     708        static int inext,inextp; 
     709        static long ma[56]; 
     710        static int iff=0; 
     711        long mj,mk; 
     712        int i,ii,k; 
     713 
     714        if (*idum < 0 || iff == 0) { 
     715                iff=1; 
     716                mj=MSEED-(*idum < 0 ? -*idum : *idum); 
     717                mj %= MBIG; 
     718                ma[55]=mj; 
     719                mk=1; 
     720                for (i=1;i<=54;i++) { 
     721                        ii=(21*i) % 55; 
     722                        ma[ii]=mk; 
     723                        mk=mj-mk; 
     724                        if (mk < MZ) mk += MBIG; 
     725                        mj=ma[ii]; 
     726                } 
     727                for (k=1;k<=4;k++) 
     728                        for (i=1;i<=55;i++) { 
     729                                ma[i] -= ma[1+(i+30) % 55]; 
     730                                if (ma[i] < MZ) ma[i] += MBIG; 
     731                        } 
     732                inext=0; 
     733                inextp=31; 
     734                *idum=1; 
     735        } 
     736        if (++inext == 56) inext=1; 
     737        if (++inextp == 56) inextp=1; 
     738        mj=ma[inext]-ma[inextp]; 
     739        if (mj < MZ) mj += MBIG; 
     740        ma[inext]=mj; 
     741        return mj*FAC; 
     742} 
     743#undef MBIG 
     744#undef MSEED 
     745#undef MZ 
     746#undef FAC 
     747 
    647748 
    648749// returns the interpolated point value in xx[0,n-1] that has the value x 
     
    696797                        return((long)Monte_SANSX2); 
    697798                        break; 
    698  
     799                case 2:                                         //  
     800                        return((long)DebyeSpheresX); 
     801                        break; 
     802                case 3:                                         //  
     803                        return((long)Monte_SANSX3); 
     804                        break; 
     805                case 4:                                         //  
     806                        return((long)Monte_SANSX4); 
     807                        break; 
    699808        } 
    700809        return(NIL); 
Note: See TracChangeset for help on using the changeset viewer.