0001 function [sampleStruct] = ACHRSamplerParallelGeneral(sampleStruct,nLoops,stepsPerPoint, maxtime, proc, fdirectory)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024 warning off MATLAB:divideByZero;
0025
0026
0027
0028
0029 if nargin < 5
0030 parallel = 0;
0031 proc = 0;
0032 elseif nargin >= 5 && proc == 1
0033 parallel = 0;
0034 proc = 0;
0035 else
0036 parallel = 1;
0037 if proc < 0
0038 proc = -proc;
0039 else
0040 numproc = proc;
0041 proc = 0;
0042
0043 delete('xxMasterfil*.mat');
0044 delete('xxRound*.mat');
0045 delete('xxRoundDonePrint*.mat');
0046 delete('xxRoundAck*.mat');
0047 delete('xxDoneRound*.mat');
0048 delete('xxDoneP*.mat');
0049 delete('xxGlobalDone*.mat*');
0050 end
0051 end
0052
0053
0054 maxMinTol = 1e-10;
0055
0056 uTol = 1e-10;
0057 safetycheck = false;
0058
0059 totalStepCount = 0;
0060 t0 = clock;
0061
0062
0063 if proc == 0
0064 if( ~ isfield(sampleStruct, 'points'))
0065 points = sampleStruct.warmupPts;
0066 else
0067 points = sampleStruct.points;
0068 end
0069 offset = sampleStruct.internal.offset;
0070
0071 [dimX,nPoints] = size(points);
0072
0073 points = points - offset*ones(1, nPoints);
0074
0075 ub = sampleStruct.internal.ubnew;
0076 lb = sampleStruct.internal.lbnew;
0077 A = sampleStruct.internal.Anew;
0078 C = sampleStruct.internal.Cnew;
0079 D = sampleStruct.internal.Dnew;
0080 fixed = union(sampleStruct.internal.fixed,find(ub==lb));
0081 if (~isfield(sampleStruct.internal,'N'))
0082 if size(A,1)==0
0083 N=[];
0084 sampleStruct.internal.N=N;
0085 else
0086 if issparse(A)
0087 N = null(full(A));
0088 else
0089 N = null(A);
0090 end
0091 sampleStruct.internal.N = N;
0092 end
0093 else
0094 N = sampleStruct.internal.N;
0095 end
0096
0097 movable = (1:dimX)';
0098 movable(fixed) = [];
0099 if safetycheck
0100 Nsmall = null(full(A(:, movable)));
0101 else
0102 Nsmall = [];
0103 end
0104
0105
0106 centerPoint = mean(points, 2);
0107
0108 fidErr = fopen('ACHRParallelError.txt','w');
0109
0110 pointRange = 1:nPoints;
0111 totalloops = nLoops;
0112
0113 if parallel
0114 blah = 1;
0115 display('saving master file.');
0116 save('xxMasterfile', 'ub', 'lb', 'A', 'C', 'D', 'fixed', 'N', 'movable', 'Nsmall', 'numproc', 'nPoints' );
0117 display('finished saving master file. spawning processes');
0118 for i = 1:(numproc - 1)
0119 command = strcat('matlab -singleCompThread -automation -nojvm -r ACHRSamplerParallelGeneral([],',num2str(nLoops),',',num2str(stepsPerPoint),',0,', num2str(-i) ,',''' ,pwd, ''');exit; &' );
0120 display(command)
0121 system(command);
0122 end
0123 display('finished spawning processes');
0124 end
0125 end
0126
0127 if proc > 0
0128 cd (fdirectory);
0129 load('xxMasterfile', 'ub', 'lb', 'A', 'C', 'D', 'fixed', 'N', 'movable', 'Nsmall', 'numproc', 'nPoints');
0130 blah = 1;
0131 end
0132
0133 for i = 1:nLoops
0134 if parallel
0135 if proc == 0
0136
0137 display(strcat('distributing points round ', num2str(i)));
0138 save(strcat('xxRound', num2str(i)), 'points', 'centerPoint');
0139 save(strcat('xxRoundDonePrint', num2str(i)), 'blah');
0140 display(strcat('finished distributing points round ', num2str(i)));
0141 else
0142 display(strcat('reading in points round ', num2str(i)));
0143 while exist(strcat('xxRoundDonePrint', num2str(i), '.mat'), 'file') ~= 2;
0144
0145 fprintf(1, '.');
0146 if exist('xxGlobalDone.mat', 'file') == 2
0147 exit;
0148 end
0149 pause(.25);
0150 end
0151 fprintf(1,'\nloading files four next round.\n');
0152 try load(strcat('xxRound', num2str(i)), 'points', 'centerPoint');
0153 catch pause(15)
0154 load(strcat('xxRound', num2str(i)), 'points', 'centerPoint');
0155 end
0156 save(strcat('xxRoundAck', num2str(i),'x', num2str(proc) ), 'blah');
0157 display(strcat('finished reading input and acknowledgment sent ', num2str(i)));
0158 end
0159
0160 pointRange = subparts(nPoints, numproc, proc);
0161 end
0162
0163
0164 for pointCount = pointRange
0165
0166 randVector = rand(stepsPerPoint,1);
0167 prevPoint = points(:,pointCount);
0168 curPoint = prevPoint;
0169 if mod(pointCount,200) == 0
0170 display(pointCount);
0171 end
0172 saveCoords = prevPoint(fixed);
0173 for stepCount = 1:stepsPerPoint
0174
0175 randPointID = ceil(nPoints*rand);
0176 randPoint = points(:,randPointID);
0177
0178
0179 u = (randPoint-centerPoint);
0180 if ~isempty(fixed)
0181
0182 if safetycheck
0183 u(movable) = Nsmall * (Nsmall' * u(movable));
0184 end
0185
0186
0187 u(fixed) = 0;
0188 end
0189 u = u/norm(u);
0190
0191
0192 distUb = (ub - prevPoint);
0193 distLb = (prevPoint - lb);
0194 distD = (D-C*prevPoint);
0195
0196
0197 posDirn = (u > uTol);
0198 negDirn = (u < -uTol);
0199 move = C*u;
0200 posDirn2 = (move > uTol);
0201 negDirn2 = (move < -uTol);
0202
0203
0204 maxStepTemp = distUb./u;
0205 minStepTemp = -distLb./u;
0206 StepD = distD./move;
0207 maxStepVec = [maxStepTemp(posDirn);minStepTemp(negDirn);StepD(posDirn2 )];
0208 minStepVec = [minStepTemp(posDirn);maxStepTemp(negDirn);StepD(negDirn2 )];
0209
0210
0211 maxStep = min(maxStepVec);
0212 minStep = max(minStepVec);
0213
0214
0215 if (abs(minStep) < maxMinTol && abs(maxStep) < maxMinTol) || (minStep > maxStep)
0216 fprintf('Warning small step: %f %f\n',minStep,maxStep);
0217 continue;
0218 end
0219
0220
0221
0222 stepDist = minStep + randVector(stepCount)*(maxStep-minStep);
0223
0224
0225
0226 curPoint = prevPoint + stepDist*u;
0227
0228
0229 if mod (stepCount, 25) == 0
0230 if ~isempty(N)
0231 curPoint = N* (N' * curPoint);
0232 end
0233 curPoint(fixed) = saveCoords;
0234 end
0235
0236
0237 if (mod(totalStepCount,1000)==0) && proc == 0
0238 fprintf(fidErr,'%10.8f\t%10.8f\t',max(curPoint-ub),max(lb-curPoint));
0239 end
0240
0241
0242 overInd = (curPoint > ub);
0243 underInd = (curPoint < lb);
0244 if (sum(overInd)>0) || (sum(underInd)>0)
0245 curPoint(overInd) = ub(overInd);
0246 curPoint(underInd) = lb(underInd);
0247 end
0248
0249
0250 if (mod(totalStepCount,1000) == 0) && proc == 0
0251 fprintf(fidErr,'%10.8f\n',full(max(max(abs(A*curPoint)))));
0252 end
0253
0254 prevPoint = curPoint;
0255
0256
0257 totalStepCount = totalStepCount + 1;
0258
0259 end
0260
0261
0262 if ~isempty(N)
0263 curPoint = N* (N' * curPoint);
0264 end
0265 curPoint(fixed) = saveCoords;
0266 centerPoint = centerPoint + (curPoint - points(:,pointCount))/nPoints;
0267
0268
0269 points(:,pointCount) = curPoint;
0270 end
0271
0272 if parallel
0273 if proc == 0
0274
0275 display(strcat ('waiting for acknowledgement', num2str(i)));
0276 donewaiting = 0;
0277 while ~donewaiting
0278 donewaiting = 1;
0279 for k = 1:(numproc-1);
0280 if exist(strcat('xxRoundAck',num2str(i), 'x', num2str(k),'.mat'), 'file') ~= 2
0281 donewaiting = 0;
0282 end
0283 end
0284
0285 fprintf(1, '.');
0286 pause(.25)
0287 end
0288
0289 fprintf(1, '\n');
0290 for k = 1:(numproc-1)
0291 delete (strcat('xxRoundAck', num2str(i), 'x', num2str(k), '.mat'));
0292 end
0293 delete(strcat('xxRound', num2str(i), '.mat'));
0294 delete(strcat('xxRoundDonePrint', num2str(i), '.mat'));
0295
0296
0297 display(strcat ('waiting for return files ', num2str(i)));
0298 donewaiting = 0;
0299 while ~donewaiting
0300 donewaiting = 1;
0301
0302 for k = 1:(numproc-1);
0303 if exist(strcat('xxDoneP',num2str(i), 'x', num2str(k),'.mat'), 'file') ~= 2
0304 donewaiting = 0;
0305 end
0306 end
0307 pause(.25)
0308 fprintf(1,'.');
0309 end
0310 fprintf(1, '\nAll processes finished. reading\n');
0311 for k = 1:(numproc-1)
0312 load (strcat('xxDoneRound', num2str(i), 'x', num2str(k)), 'points2');
0313 r2 = subparts(nPoints, numproc, k);
0314 points(:,r2) = points2;
0315 delete (strcat('xxDoneRound', num2str(i), 'x', num2str(k), '.mat') );
0316 delete (strcat('xxDoneP',num2str(i), 'x', num2str(k), '.mat'));
0317 end
0318 centerPoint = mean(points, 2);
0319 display(strcat('done with round ', num2str(i)));
0320 else
0321 points2 = points(:, pointRange);
0322 save (strcat('xxDoneRound', num2str(i), 'x', num2str(proc)), 'points2');
0323 save (strcat('xxDoneP',num2str(i), 'x', num2str(proc)), 'blah');
0324 end
0325 end
0326
0327 t1 = clock();
0328 fprintf('%10.0f s %d steps\n',etime(t1, t0),i*stepsPerPoint);
0329 if etime(t1, t0) > maxtime && proc == 0
0330 totalloops = i;
0331 break;
0332 end
0333 end
0334
0335 if proc > 0
0336 return;
0337 end
0338 points = points + offset*ones(1, nPoints);
0339 sampleStruct.points = points;
0340 if ~ isfield(sampleStruct, 'steps')
0341 sampleStruct.steps = 0;
0342 end
0343
0344 sampleStruct.steps = sampleStruct.steps + stepsPerPoint*totalloops;
0345
0346 if parallel
0347 save('xxGlobalDone', 'blah');
0348 delete('xxMaster*.mat');
0349 end
0350
0351 fclose(fidErr);
0352
0353 function out = subparts(nPoints, n, k)
0354 out = (floor(nPoints*k/n)+1) : (floor(nPoints*(k+1)/n));
0355 return