Revolving Doors Riddle - Matlab Time-Efficient Sparse Matrix Use

I'm running a code with many iterations using large sparse matrices. There are three lines in my code that take about 75% of the running time and I think I can use the special structure of my sparse matrix to reduce that time, but so far I haven't managed to do it. I would love your help!!

Ok, here's the gist of my code:

I = 70;
J = 1000;
A = rand(I);
A = A./repmat(sum(A, 2), 1, I);
S   = kron(A, speye(J));
indj = randi(J,I,1); 
tic
for i = 1:I
    S(:, (i-1)*J+indj(i)) = sum(S(:, (i-1)*J + (1:indj(i))), 2);
end
toc

You can skip the following 2 paragraphs

Here's a story to make the example a bit more lively. An old man is visiting sick people at different hospitals. There are 1000 (J) hospitals, and each hospital has 70 (I) rooms in it. The matrix A is the transition matrix that specifies the probability of the old man moving from one room at the hospital to another room within the same hospital. A(i1,i2) is the probability the old man moves from room i1 to room i2 (so columns sum to 1). The big S matrix is the transition probability matrix, where moving from room i1 at hospital j1 to room i2 at hospital j2 is given by the (J*(i1-1)+j1, J*(i2-1)+j2) element. There is no way the old man moves from one hospital to another, so the matrix is sparse.

Something magical happens and now all the doors to room number i in the first indj(i) hospitals all lead to the same hospital, hospital indj(i). So the old man can now magically move between hospitals. We need to change the S matrix accordingly. This amounts to two things, increasing the probability of moving to room i at hospital indj(i), for all i, and setting to zero the probability of getting into all rooms lower than indj(i) at hospital i, for all i. The latter I can do very efficiently, but the first part is taking me too long.

Why I think there's a chance to reduce running time

1. Loop. The part between the tic and toc can be written without a loop. I have done it, but it made it run much slower perhaps because the length of the sub2ind is very large.

2. Matrix structure. Notice that we don’t need the entire sum, only one element needs to be added. These loops achieve the same outcome (but here, obviously, much slower):

   for i = 1:I
   for ii = 1:I
   for j = 1:indj(i)-1
       S((ii-1)*J+j, (i-1)*J+indj(i)) = S((ii-1)*J+j, (i-1)*J+indj(i)) + S((ii-1)*J+j, (i-1)*J+j);
   end
   end
   end
   

This makes me somewhat hopeful that there is a way to make the calculation faster…

Your help is HIGHLY appreciated!